Regression Project¶

In this project, your goal is to build regression models of housing prices. The models should learn from data and be able to predict the median house price in a district (which is a population of 600 to 3000 people), given some predictor variable

Das Vorgehen in diesem Notebook orientiert sich am Data Science Lifecycle. Es wird eine Modellpräzision von > 70% gewünscht. Das bedeutet, dass der R-squared Wert des besten Modells >0,7 sein soll.

Python Setup¶

Zunächst werden mit der Import Methode wichtige Python Bibliotheken importiert.

import pandas as pd
import numpy as np 

import seaborn as sns
sns.set_theme(style='ticks')

import matplotlib.pyplot as plt
%matplotlib inline

import statsmodels.formula.api as smf
import statsmodels.api as sm
import statsmodels.gam.api as smg

from sklearn.compose import ColumnTransformer
from sklearn.compose import make_column_selector as selector
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
from sklearn.model_selection import train_test_split
import sklearn.metrics as skm
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import Lasso
from sklearn.preprocessing import SplineTransformer
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline

import plotly.express as px

Datenaufbereitung¶

Daten laden¶

Importieren der Daten aus der CSV in Github und Erstellung eines Pandas Dataframes zur Bearbeitung und Verwendung der Daten in diesem Notebook.

GIT = "https://raw.githubusercontent.com/jan-kirenz/project-annikas428/main/"
DATA = "project_data.csv"
TOKEN = "?token=GHSAT0AAAAAABPURMBOEW2WS3JEBXGI5TPAYPJRBIQ"

df_initial = pd.read_csv(GIT + DATA + TOKEN)

Inspect Data¶

Im ersten Schritt ist es wichtig, sich einen Überblick über die Daten zu verschaffen. Begonnen wird mit einer Anzeige einiger Datensätze.

df_initial

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	price_category
0	-122.23	37.88	41.0years	880	129.0	322	126	8.3252	452600.0$	NEAR BAY	above
1	-122.22	37.86	21.0	7099	1106.0	2401	1138	8.3014	358500.0	NEAR BAY	above
2	-122.24	37.85	52.0	1467	190.0	496	177	7.2574	352100.0	NEAR BAY	above
3	-122.25	37.85	52.0	1274	235.0	558	219	5.6431	341300.0	NEAR BAY	above
4	-122.25	37.85	52.0	1627	280.0	565	259	3.8462	342200.0	NEAR BAY	above
...	...	...	...	...	...	...	...	...	...	...	...
20635	-121.09	39.48	25.0	1665	374.0	845	330	1.5603	78100.0	INLAND	above
20636	-121.21	39.49	18.0	697	150.0	356	114	2.5568	77100.0	INLAND	above
20637	-121.22	39.43	17.0	2254	485.0	1007	433	1.7000	92300.0	INLAND	above
20638	-121.32	39.43	18.0	1860	409.0	741	349	1.8672	84700.0	INLAND	above
20639	-121.24	39.37	16.0	2785	616.0	1387	530	2.3886	89400.0	INLAND	above

20640 rows × 11 columns

Beim Sichten der Datensätze fällt auf, dass im ersten Datensatz zwei Unstimmigkeiten im Vergleich zu den anderen Datensätzen vorliegen:
median_house_value ist nur hier mit der Währung Dollar angegeben und hinter housing_median_age steht der Zusatz “years”.

Nun folgt Pandas Dataframe Info Methode. Diese zeigt die Variablen und deren Datentypen, sowie die Anzahl der non-null values je Spalte an.

df_initial.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 20640 entries, 0 to 20639
Data columns (total 11 columns):
 #   Column              Non-Null Count  Dtype  
---  ------              --------------  -----  
 0   longitude           20640 non-null  float64
 1   latitude            20640 non-null  float64
 2   housing_median_age  20640 non-null  object 
 3   total_rooms         20640 non-null  int64  
 4   total_bedrooms      20433 non-null  float64
 5   population          20640 non-null  int64  
 6   households          20640 non-null  int64  
 7   median_income       20640 non-null  float64
 8   median_house_value  20640 non-null  object 
 9   ocean_proximity     20640 non-null  object 
 10  price_category      20640 non-null  object 
dtypes: float64(4), int64(3), object(4)
memory usage: 1.7+ MB

Die Variablen housing_median_age und und median_house_value wurden als Datentyp object definiert. Dies muss genauer untersucht werden, da bei einer Jahreszahl bzw. einem monetären Wert ein numerischer Datentyp (int oder float) angemessener ist. Die Ursache hierfür ist auf die Zusatzangaben “$” und “years” zurückzuführen.
Die Varialben ocean_proximity und price_category sind keine numerischen sondern kategoriale Datentypen. Diese lassen sich jeweils in einer internen Rangordnung anordnen. Sie sind daher ordinal. Der richtige Datentyp hierfür ist category. Die Transformation der Datentypen erfolgt in Kapitel “Data Transformation”.
Neben den Datentypen wird auch die Anzahl der non-null values je Variable angezeigt. Hier haben alle Variablen eine Anzahl von 20640 non-null values. Lediglich die Variable total_bedrooms liegt unter diesem Wert. Dies lässt darauf schließen, dass hier einige Datensätze keinen Wert haben. Dies gilt es im Folgenden genauer zu untersuchen.

# show missing values -if present - in yellow
sns.heatmap(df_initial.isnull(), 
            yticklabels=False,
            cbar=False, 
            cmap='viridis');

print(df_initial.isnull().sum())

longitude               0
latitude                0
housing_median_age      0
total_rooms             0
total_bedrooms        207
population              0
households              0
median_income           0
median_house_value      0
ocean_proximity         0
price_category          0
dtype: int64

Wie bereits vermutet, enthält die Variable total_bedrooms 207 null values. Die Abbildung visualisiert diese nicht befüllten Datensätze und deren Verteilung im Datenset in gelb. Die null-values sind zufällig im Datenset verteilt. Da das Datenset über 20000 Datensätze enthält, werden die Zeilen mit den null-values im Kapitel Transform Data entfernt. Alternativ wäre auch eine Belegung der null-values beispielsweise mit dem Durchschnittswert der Variable möglich. Darauf wird hier jedoch verzichtet, genug andere Datensätze vorliegen.

Datentransformation¶

Entfernung der null values¶

# drop rows with missing values
df = df_initial.dropna()
print(df.isnull().sum())

longitude             0
latitude              0
housing_median_age    0
total_rooms           0
total_bedrooms        0
population            0
households            0
median_income         0
median_house_value    0
ocean_proximity       0
price_category        0
dtype: int64

Entfernung unwichtiger Features¶

Das Feature -price_category_ wird entfernt, da es sich direkt aus der Variablen median_house_value ableitet und daher nicht zur Berechnung des Wertes genutzt werden darf. Dieses Feature wird im Notebook Classification als abhängige Variable genutzt.

#drop irrelevant features
df = df.drop(['price_category'], axis=1)

Konvertierung von Datentypen¶

Die mithilfe der info() Methode identifizierten Features, welche eine Datentyp Konvertierung benötigen, werden nun bearbeitet.

# Convert data types Ocean_Proximity and Price_Category to Data Type "Category" as they are categorical variables (ordinal)
df.ocean_proximity = df.ocean_proximity.astype('category')

# Housing_Median_Age and Median_House_Value are actually int64 values, but pandas did not get this because of signs like '€'
df.housing_median_age = df.housing_median_age.str.replace('years', '', regex=False)
df.median_house_value = df.median_house_value.str.replace('$', '', regex=False)
df.housing_median_age = pd.to_numeric(df.housing_median_age, errors='coerce').astype('int64')
df.median_house_value = pd.to_numeric(df.median_house_value, errors='coerce').astype('float64')
df.total_bedrooms = df.total_bedrooms.astype('int64')
df.info()
df.head()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 20433 entries, 0 to 20639
Data columns (total 10 columns):
 #   Column              Non-Null Count  Dtype   
---  ------              --------------  -----   
 0   longitude           20433 non-null  float64 
 1   latitude            20433 non-null  float64 
 2   housing_median_age  20433 non-null  int64   
 3   total_rooms         20433 non-null  int64   
 4   total_bedrooms      20433 non-null  int64   
 5   population          20433 non-null  int64   
 6   households          20433 non-null  int64   
 7   median_income       20433 non-null  float64 
 8   median_house_value  20433 non-null  float64 
 9   ocean_proximity     20433 non-null  category
dtypes: category(1), float64(4), int64(5)
memory usage: 1.6 MB

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity
0	-122.23	37.88	41	880	129	322	126	8.3252	452600.0	NEAR BAY
1	-122.22	37.86	21	7099	1106	2401	1138	8.3014	358500.0	NEAR BAY
2	-122.24	37.85	52	1467	190	496	177	7.2574	352100.0	NEAR BAY
3	-122.25	37.85	52	1274	235	558	219	5.6431	341300.0	NEAR BAY
4	-122.25	37.85	52	1627	280	565	259	3.8462	342200.0	NEAR BAY

Mithilfe der string replace Methode werrden die zuvor identifizierten Zeichen “years” und “$” entfernt. Anschließend ist eine Transformation der Features in integer bzw. float Werte möglich. Für housing_median_age wird integer gewählt, da die Jahreszahlen als Ganze Zahlen also discrete Werte angegeben sind.
Für median_house_value wird der Datentyp float benötigt, da die Ergebnisse einer linearen Regression stehts continous sind.

Hinzufügen neuer Features¶

df = df.assign(people_per_household=lambda df: df.population / df.households)
df = df.assign(rooms_per_household=lambda df: df.total_rooms / df.households)
df = df.assign(bedrooms_per_people=lambda df: df.total_bedrooms / df.population)
df = df.assign(bedrooms_per_rooms=lambda df: df.total_bedrooms / df.total_rooms)

Diese Features werden dem Datenset hinzugefügt, in der Hoffnung, dass diese für die Verwendung in den Modellen nützlich sein könnten. Untersucht wird dies im Kapitel EDA.

Datensplit¶

Split des Datensets in 80% Trainingsdaten und 20% Testdaten. Die Testdaten werden anschließend bis zur finalen Evaluation der besten Modelle nicht verwendet.

train_dataset = df.sample(frac=0.8, random_state=0) 
test_dataset = df.drop(train_dataset.index)

train_dataset

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
14185	-117.08	32.70	37	2176	418	1301	375	2.8750	98900.0	NEAR OCEAN	3.469333	5.802667	0.321291	0.192096
6125	-117.91	34.11	20	3158	684	2396	713	3.5250	153000.0	<1H OCEAN	3.360449	4.429173	0.285476	0.216593
14095	-117.10	32.75	11	2393	726	1905	711	1.3448	91300.0	NEAR OCEAN	2.679325	3.365682	0.381102	0.303385
14359	-117.22	32.74	52	1260	202	555	209	7.2758	345200.0	NEAR OCEAN	2.655502	6.028708	0.363964	0.160317
18004	-121.99	37.29	32	2930	481	1336	481	6.4631	344100.0	<1H OCEAN	2.777547	6.091476	0.360030	0.164164
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
6168	-117.93	34.09	35	1891	353	1093	382	4.0167	165500.0	<1H OCEAN	2.861257	4.950262	0.322964	0.186674
4174	-118.21	34.10	36	2000	533	1234	535	3.7437	241700.0	<1H OCEAN	2.306542	3.738318	0.431929	0.266500
3593	-118.48	34.24	32	2621	412	1285	414	6.6537	267600.0	<1H OCEAN	3.103865	6.330918	0.320623	0.157192
5366	-118.39	34.04	44	1873	286	635	283	5.5951	461300.0	<1H OCEAN	2.243816	6.618375	0.450394	0.152696
17879	-122.00	37.40	17	5121	1017	1470	968	2.9706	81300.0	<1H OCEAN	1.518595	5.290289	0.691837	0.198594

16346 rows × 14 columns

Deskriptive Statistik¶

Kennzahlenanalyse numerischer Variablen¶

Mithilfe der describe() Methode werden die wichtigsten Kennzahlen numerischer Variablen (bspw. Durchschnitt, min und max Werte) ausgegeben.

#summary statistics for all numerical columns
train_dataset.describe().T

	count	mean	std	min	25%	50%	75%	max
longitude	16346.0	-119.564154	2.002618	-124.350000	-121.790000	-118.490000	-118.000000	-114.470000
latitude	16346.0	35.630318	2.138574	32.550000	33.930000	34.250000	37.710000	41.950000
housing_median_age	16346.0	28.664505	12.556764	1.000000	18.000000	29.000000	37.000000	52.000000
total_rooms	16346.0	2622.235776	2169.548287	11.000000	1448.000000	2119.000000	3120.750000	39320.000000
total_bedrooms	16346.0	535.281659	418.469078	3.000000	296.000000	432.500000	644.000000	6445.000000
population	16346.0	1416.087055	1103.842065	3.000000	784.250000	1164.000000	1711.000000	28566.000000
households	16346.0	496.758167	379.109535	3.000000	280.000000	408.000000	600.000000	6082.000000
median_income	16346.0	3.869337	1.902228	0.499900	2.555675	3.533200	4.744225	15.000100
median_house_value	16346.0	206916.154411	115676.394484	14999.000000	119300.000000	179700.000000	265900.000000	500001.000000
people_per_household	16346.0	3.089529	11.525259	0.692308	2.428850	2.816514	3.282584	1243.333333
rooms_per_household	16346.0	5.423747	2.160962	0.846154	4.443697	5.230149	6.051703	62.422222
bedrooms_per_people	16346.0	0.403370	0.223233	0.000670	0.315480	0.372646	0.443446	8.750000
bedrooms_per_rooms	16346.0	0.213188	0.057891	0.100000	0.175591	0.203263	0.239668	1.000000

Der Durchschnitt (mean) der Variablen housing_median_age über alle Datensätze beträgt 28,64 Jahre. Die Standardabweichung, also die durchschnittliche Abweichung der Werte vom Durchschnitt, beträgt 12,56 Jahre. Das Minimum liegt bei 1 Jahr, das Maximum der 52 Jahren. 25% der Datensätze liegen unter 18 Jahren (1. Quartil). 50% zwischen unter 29 Jahren (2. Quartil) und 75% zwischen unter 37 Jahren (3. Quartil).
Im Durchschnitt beträgt der Wert people_per_household 3,09. Die durchschnittliche Abweichung von diesem Wert beträgt 11,53. Das Minimum liegt bei 0,69 und das Maximum bei 1243 Personen pro Haushalt. 25% der Datensätze liegen zwischen 0,69 und 2,43. 50% unter 2,82 und 75% der Datensätze unter 3,28 Personen pro Haushalt.

print(train_dataset.quantile(0.99, interpolation='nearest'))

longitude                 -116.290000
latitude                    40.630000
housing_median_age          52.000000
total_rooms              11021.000000
total_bedrooms            2208.000000
population                5723.000000
households                1955.000000
median_income               10.557500
median_house_value      500001.000000
people_per_household         5.286920
rooms_per_household         10.500000
bedrooms_per_people          0.941606
bedrooms_per_rooms           0.407472
Name: 0.99, dtype: float64

Die Quantile() Methode gibt für jede Variable einen Wert an, unter welchem (in diesem Beispiel) 99% der Datensätze liegen.
Beispielsweise liegen 99% der Datensätze unter einer Anzahl total_rooms in Höhe von 11021.
99% der Distrikte haben eine population unter 5723 Einwohnern.
Im Nachfolgenden werden einige mögliche Ausreiser identifiziert.

Dieser Datensatz steht für ein Distrikt, welches aus nur 1 Haushalt besteht. In diesem Haushalt leben 13 Leute, die sich 8 Räume, wovon einer ein Schlafzimmer ist, teilen. Der Wert des Hauses beträgt 500001 $.
Dies erscheint etwas unwahrscheinlich. Womöglich handelt es sich um eine Falscheingabe. Der Datensatz wird daher entfernt.

train_dataset.loc[train_dataset.people_per_household > 100]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
19006	-121.98	38.32	45	19	5	7460	6	10.2264	137500.0	INLAND	1243.333333	3.166667	0.000670	0.263158
16669	-120.70	35.32	46	118	17	6532	13	4.2639	350000.0	NEAR OCEAN	502.461538	9.076923	0.002603	0.144068
3364	-120.51	40.41	36	36	8	4198	7	5.5179	67500.0	INLAND	599.714286	5.142857	0.001906	0.222222

Mehr als 50 Leute in einem Haushalten wirken sehr unrealistisch. 99% der Datensätze liegen unter einem Wert von 5,39 Personen pro Haushalt. Daher werden diese Datensätze als Ausreiser identifiziert.

train_dataset.loc[train_dataset.bedrooms_per_people > 8]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
11862	-121.25	40.27	25	958	245	28	16	2.625	67500.0	INLAND	1.75	59.875	8.75	0.255741

Auch diese Werte liegen weit über dem 99% Perzentil der Varialben bedrooms_per_people und werden daher als Ausreiser entfernt. Zwei dieser Werte sind auch bereits rooms_per_household Ausreiser.

train_dataset.loc[train_dataset.rooms_per_household > 55]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
12447	-114.49	33.97	17	2809	635	83	45	1.6154	87500.0	INLAND	1.844444	62.422222	7.650602	0.226059
11862	-121.25	40.27	25	958	245	28	16	2.6250	67500.0	INLAND	1.750000	59.875000	8.750000	0.255741
1913	-120.06	39.01	19	2967	528	112	48	4.0714	437500.0	INLAND	2.333333	61.812500	4.714286	0.177958

#drop people per household outliers
train_dataset = train_dataset.drop([3364, 16669, 19006])
#drop bedroooms per people outliers
train_dataset = train_dataset.drop([11862])
#drop rooms per household_outliers
train_dataset = train_dataset.drop([1913, 12447])

#summary statistics for all numerical columns
train_dataset.describe().T

	count	mean	std	min	25%	50%	75%	max
longitude	16340.0	-119.564056	2.002423	-124.350000	-121.790000	-118.490000	-118.000000	-114.470000
latitude	16340.0	35.629491	2.138024	32.550000	33.930000	34.250000	37.710000	41.950000
housing_median_age	16340.0	28.663525	12.556964	1.000000	18.000000	29.000000	37.000000	52.000000
total_rooms	16340.0	2622.775949	2169.626947	11.000000	1449.000000	2119.500000	3121.250000	39320.000000
total_bedrooms	16340.0	535.390208	418.478468	3.000000	296.000000	433.000000	644.000000	6445.000000
population	16340.0	1415.480171	1101.940296	3.000000	785.000000	1164.000000	1711.000000	28566.000000
households	16340.0	496.932313	379.070022	3.000000	280.000000	408.000000	600.250000	6082.000000
median_income	16340.0	3.869024	1.901773	0.499900	2.555600	3.533000	4.744075	15.000100
median_house_value	16340.0	206921.906977	115662.820190	14999.000000	119300.000000	179700.000000	265900.000000	500001.000000
people_per_household	16340.0	2.946757	1.291597	0.692308	2.428869	2.816514	3.282433	83.171429
rooms_per_household	16340.0	5.413407	2.023690	0.846154	4.443662	5.229614	6.050922	52.848214
bedrooms_per_people	16340.0	0.402226	0.202989	0.011222	0.315537	0.372646	0.443428	7.925926
bedrooms_per_rooms	16340.0	0.213188	0.057896	0.100000	0.175588	0.203260	0.239666	1.000000

Analyse kategorialer Variablen¶

# summary statistics for all categorical columns
train_dataset.describe(include=['category']).transpose()

	count	unique	top	freq
ocean_proximity	16340	5	<1H OCEAN	7216

Die Variable ocean_proximity besitzt 5 Ausprägungen, wovon die Ausprägung <1H OCEAN mit einer Anzhl von 7215 am häufigsten vorkommt.

train_dataset['ocean_proximity'].value_counts()

<1H OCEAN     7216
INLAND        5216
NEAR OCEAN    2107
NEAR BAY      1796
ISLAND           5
Name: ocean_proximity, dtype: int64

Hier wird die Häufigkeit einzelner Ausprägungen von ocean_proximity dargestellt. Die Ausprägung ISLAND kommt über das gesamte Datenset mit über 20000 Datensätzen nur fünfmal vor. In Relation zur Gesamtanzahl der Datensätzen ist dies sehr wenig, wodurch die Aussagekraft dieser Ausprägung vermutlich eher gering ist. Daher werden die Datensätze mit Ausprägung ISLAND mit der Ausprägung NEAR OCEAN vereint.

train_dataset.loc[train_dataset.ocean_proximity == 'ISLAND']

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
8318	-118.48	33.43	29	716	214	422	173	2.6042	287500.0	ISLAND	2.439306	4.138728	0.507109	0.298883
8316	-118.32	33.33	52	2127	512	733	288	3.3906	300000.0	ISLAND	2.545139	7.385417	0.698499	0.240715
8315	-118.33	33.34	52	2359	591	1100	431	2.8333	414700.0	ISLAND	2.552204	5.473318	0.537273	0.250530
8314	-118.32	33.35	27	1675	521	744	331	2.1579	450000.0	ISLAND	2.247734	5.060423	0.700269	0.311045
8317	-118.32	33.34	52	996	264	341	160	2.7361	450000.0	ISLAND	2.131250	6.225000	0.774194	0.265060

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR OCEAN',
    'NEAR BAY': 'NEAR BAY',
    'ISLAND': 'NEAR OCEAN',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Exploratory Data Anaysis¶

Dieses Kapitel dient der ausführlichen Datenanalyse der Testdaten. Hierfür werden vor allem geeignete Visualisierungen in Form von Diagrammen verwendet.

erster Überblick¶

# summary statistics for all numerical columns
round(train_dataset.describe(),2).transpose()

	count	mean	std	min	25%	50%	75%	max
longitude	16340.0	-119.56	2.00	-124.35	-121.79	-118.49	-118.00	-114.47
latitude	16340.0	35.63	2.14	32.55	33.93	34.25	37.71	41.95
housing_median_age	16340.0	28.66	12.56	1.00	18.00	29.00	37.00	52.00
total_rooms	16340.0	2622.78	2169.63	11.00	1449.00	2119.50	3121.25	39320.00
total_bedrooms	16340.0	535.39	418.48	3.00	296.00	433.00	644.00	6445.00
population	16340.0	1415.48	1101.94	3.00	785.00	1164.00	1711.00	28566.00
households	16340.0	496.93	379.07	3.00	280.00	408.00	600.25	6082.00
median_income	16340.0	3.87	1.90	0.50	2.56	3.53	4.74	15.00
median_house_value	16340.0	206921.91	115662.82	14999.00	119300.00	179700.00	265900.00	500001.00
people_per_household	16340.0	2.95	1.29	0.69	2.43	2.82	3.28	83.17
rooms_per_household	16340.0	5.41	2.02	0.85	4.44	5.23	6.05	52.85
bedrooms_per_people	16340.0	0.40	0.20	0.01	0.32	0.37	0.44	7.93
bedrooms_per_rooms	16340.0	0.21	0.06	0.10	0.18	0.20	0.24	1.00

import plotly.graph_objects as go

fig = go.Figure(data=go.Scattergeo(
        lon = train_dataset['longitude'],
        lat = train_dataset['latitude'],
        text = 'median house value in district: ' + train_dataset['median_house_value'].astype(str) + '$',
        marker_color = train_dataset['median_house_value'], marker_size=2, marker_colorscale = 'Inferno', marker_colorbar= {'title': 'Median House Value'}
        ))

fig.update_layout(
        title = 'Median House Value per District in California',
        geo_scope='usa',
    )
fig.show()

Die Datensätze stammen alle aus Kalifornien, USA.
Die Werte Longitude und Latitude unterscheiden sich daher kaum, wie sich auch bereits im Kapitel deskriptive Statistik gezeigt hat. Daher sind sie als Feature eher ungeeignet.
Allerdings scheint die Entfernung zum Ozean eine wichtige Rolle im Bezug auf den median_house_value zu spielen. Ocean_proximity scheint daher als feature interessant und sollte genauer untersucht werden.

train_dataset = train_dataset.drop(['longitude', 'latitude'], axis = 1)

sns.pairplot(train_dataset[['housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_house_value']])

<seaborn.axisgrid.PairGrid at 0x1ec11b17490>

Error in callback <function flush_figures at 0x000001EC10073040> (for post_execute):

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
C:\ProgramData\Anaconda3\lib\site-packages\matplotlib_inline\backend_inline.py in flush_figures()
    119         # ignore the tracking, just draw and close all figures
    120         try:
--> 121             return show(True)
    122         except Exception as e:
    123             # safely show traceback if in IPython, else raise

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib_inline\backend_inline.py in show(close, block)
     39     try:
     40         for figure_manager in Gcf.get_all_fig_managers():
---> 41             display(
     42                 figure_manager.canvas.figure,
     43                 metadata=_fetch_figure_metadata(figure_manager.canvas.figure)

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\display.py in display(include, exclude, metadata, transient, display_id, *objs, **kwargs)
    318             publish_display_data(data=obj, metadata=metadata, **kwargs)
    319         else:
--> 320             format_dict, md_dict = format(obj, include=include, exclude=exclude)
    321             if not format_dict:
    322                 # nothing to display (e.g. _ipython_display_ took over)

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\formatters.py in format(self, obj, include, exclude)
    178             md = None
    179             try:
--> 180                 data = formatter(obj)
    181             except:
    182                 # FIXME: log the exception

C:\ProgramData\Anaconda3\lib\site-packages\decorator.py in fun(*args, **kw)
    230             if not kwsyntax:
    231                 args, kw = fix(args, kw, sig)
--> 232             return caller(func, *(extras + args), **kw)
    233     fun.__name__ = func.__name__
    234     fun.__doc__ = func.__doc__

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\formatters.py in catch_format_error(method, self, *args, **kwargs)
    222     """show traceback on failed format call"""
    223     try:
--> 224         r = method(self, *args, **kwargs)
    225     except NotImplementedError:
    226         # don't warn on NotImplementedErrors

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

C:\ProgramData\Anaconda3\lib\site-packages\IPython\core\pylabtools.py in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
    149         FigureCanvasBase(fig)
    150 
--> 151     fig.canvas.print_figure(bytes_io, **kw)
    152     data = bytes_io.getvalue()
    153     if fmt == 'svg':

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
   2292             if bbox_inches:
   2293                 if bbox_inches == "tight":
-> 2294                     bbox_inches = self.figure.get_tightbbox(
   2295                         renderer, bbox_extra_artists=bbox_extra_artists)
   2296                     if pad_inches is None:

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\figure.py in get_tightbbox(self, renderer, bbox_extra_artists)
   1626 
   1627         for a in artists:
-> 1628             bbox = a.get_tightbbox(renderer)
   1629             if bbox is not None and (bbox.width != 0 or bbox.height != 0):
   1630                 bb.append(bbox)

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\artist.py in get_tightbbox(self, renderer)
    353             The enclosing bounding box (in figure pixel coordinates).
    354         """
--> 355         bbox = self.get_window_extent(renderer)
    356         if self.get_clip_on():
    357             clip_box = self.get_clip_box()

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\patches.py in get_window_extent(self, renderer)
    621 
    622     def get_window_extent(self, renderer=None):
--> 623         return self.get_path().get_extents(self.get_transform())
    624 
    625     def _convert_xy_units(self, xy):

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\path.py in get_extents(self, transform, **kwargs)
    612         from .transforms import Bbox
    613         if transform is not None:
--> 614             self = transform.transform_path(self)
    615         if self.codes is None:
    616             xys = self.vertices

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in transform_path(self, path)
   1606         that began as line segments.
   1607         """
-> 1608         return self.transform_path_affine(self.transform_path_non_affine(path))
   1609 
   1610     def transform_path_affine(self, path):

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in transform_path_affine(self, path)
   1616         ``transform_path_affine(transform_path_non_affine(values))``.
   1617         """
-> 1618         return self.get_affine().transform_path_affine(path)
   1619 
   1620     def transform_path_non_affine(self, path):

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in get_affine(self)
   2444             return self._b.get_affine()
   2445         else:
-> 2446             return Affine2D(np.dot(self._b.get_affine().get_matrix(),
   2447                                    self._a.get_affine().get_matrix()))
   2448 

C:\ProgramData\Anaconda3\lib\site-packages\matplotlib\transforms.py in __init__(self, matrix, **kwargs)
   1912             # A bit faster than np.identity(3).
   1913             matrix = IdentityTransform._mtx.copy()
-> 1914         self._mtx = matrix.copy()
   1915         self._invalid = 0
   1916 

KeyboardInterrupt: 

Die hier abgebildeten, unabhängigen Variablen weisen eine starke Korrelation untereinander auf (bis auf housing_median_age).
In Bezug auf die abhängige Variable median_house_value zeigen die Plots zunächst keinen starken Zusammenhang.
Interessant sieht das Histogramme von median_house_value aus, da es am äußeren rechten Rand nochmal hohe Häufigkeiten haben. Dieses Histogramme wird später nochmal genauer analysisiert.

sns.pairplot(train_dataset[['median_income', 'people_per_household', 'rooms_per_household', 'bedrooms_per_people', 'bedrooms_per_rooms', 'median_house_value']]);

Zwischen median_income und median_house_value ist ein Zusammenhang erkennbar. Auch das Histogramm von median_income zeigt eine hohe Säule am rechten Rand.
Auch für die selbst definierten Variablen ist es schwierig, ein linearen Zusammenhang zu median_house_value zu erkennen.

fig = sns.displot(data= train_dataset, x='median_house_value')

Die Verteilung von Median_House_value nicht normalverteilt, rechts im hohen Preissegment ist die Häufigkeit nochmal extrem hoch. Über 500000$ exisiteren keine weiteren Hauswerte. Möglicherweise wurde der Datensatz mit Werten über 500000$ abgeschnitten und für alle höheren Werte, der Wert 500000$ gesetzt.

fig = sns.displot(data= train_dataset, x='median_income')

Womöglich wurde auch die Variable median_income nach oben abgeschnitten.

Korrelationsanalyse¶

# Create correlation matrix for numerical variables
corr_matrix = train_dataset.corr()
corr_matrix['median_house_value'].sort_values(ascending=False)

median_house_value      1.000000
median_income           0.692088
rooms_per_household     0.176506
total_rooms             0.130850
housing_median_age      0.105574
bedrooms_per_people     0.077531
households              0.063122
total_bedrooms          0.047633
population             -0.026782
longitude              -0.043011
latitude               -0.147308
people_per_household   -0.157646
bedrooms_per_rooms     -0.253006
Name: median_house_value, dtype: float64

# Make a pretty heatmap

# Use a mask to plot only part of a matrix
mask = np.zeros_like(corr_matrix)
mask[np.triu_indices_from(mask)]= True

# Change size
plt.subplots(figsize=(11, 15))

# Build heatmap with additional options
heatmap = sns.heatmap(corr_matrix, 
                      mask = mask, 
                      square = True, 
                      linewidths = .5,
                      cmap = 'coolwarm',
                      cbar_kws = {'shrink': .6,
                                'ticks' : [-1, -.5, 0, 0.5, 1]},
                      vmin = -1,
                      vmax = 1,
                      annot = True,
                      annot_kws = {"size": 10})

Laut der Korrelationsmatrix einigen sich vor allem folgende Variablen:
median_income weist mit einem Wert von 0,69 die höchste Korrelation auf, danach folgen bedrooms_per_rooms mit einer negativen Korrelation in Höhe von -0,25 und anschließend people_per_household ebenfalls mit einer negativen Korrelation in Höhe von -0,16.
Danach folgt rooms_per_household, allerdings kann diese Variable nicht in Kombination mit den beiden zuvor genannten genutzt werden, da sie sich aus total_rooms und households zusammensetzt. Diese Werte werden jedoch bereits in bedrooms_per_rooms und people_per_household verwendet, daher ist die zusätzliche Verwendung von rooms_per_household nicht sinnvoll.
housing_median_age könnte mit einer Korrelation von 0,11 noch einen positiven Einfluss auf das Modell haben.

Analyse der kategorialen Variable Ocean_Proximity¶

sns.boxenplot(data=train_dataset, x='ocean_proximity', y='median_house_value')

<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>

sns.stripplot(data=train_dataset, x='ocean_proximity', y='median_house_value', size=2)

<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>

Die Ausprägungen NEAR BAY und NEAR OCEAN unterscheiden sich hinsichtlich median_house_value kaum. Daher können diese beiden Ausprägungen zusammengefasst werden.

sns.displot(data=train_dataset, x='median_house_value', hue='ocean_proximity', kind='kde')

<seaborn.axisgrid.FacetGrid at 0x230d8270f40>

Auch in der Häufigkeitsverteilung zeigen NEAR BAY und NEAR OCEAN kaum Unterschiede.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': '<1H OCEAN' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

sns.boxplot(data=train_dataset, x='ocean_proximity', y='median_house_value')

<AxesSubplot:xlabel='ocean_proximity', ylabel='median_house_value'>

Auch NEAR COAST und <1H OCEAN zeigen eine ähnliche Verteilung des Medians und der Quartile. Das 3. Quartil von NEAR COAST zieht sich jedoch noch um ca. 50000$ weiter nach oben, als das von <1H OCEAN, daher bleibt die Differenzierung dieser beiden Ausprägungen noch bestehen.
Bei der Ausprägung INLAND sind viele Ausreiser nach oben zu erkennen.

Analyse numerischer Variablen¶

# check relationship with a joint plot
sns.jointplot(x="median_income", y="median_house_value", hue="ocean_proximity", data=train_dataset);

# check relationship with a joint plot
sns.lmplot(x="median_income", y="median_house_value", data=train_dataset);

Wie bereits aus der Korrelationsmatrix zu schließen, ist ein linearer Zusammenhang zwischen median_income und median_house_value erkennbar. Die “Linie” mit Datenpunkten bei 500000$ ergibt sich aus der anzunehmenden Begrenzung der median_house_value_ auf 500000$.
Auch innerhalb der einzelenden Ausprägungen von ocean_proximity ist eine lineare Korrelation von median_income und median_house_value erkennbar.

sns.relplot(data=train_dataset, x='total_rooms', y='median_house_value', hue='ocean_proximity', col='ocean_proximity')

<seaborn.axisgrid.FacetGrid at 0x193868f47c0>

Ein linearer Zusammenhang von total_rooms zu median_house_value ist nicht wirklich zu erkennen. Die Variable ist daher für eine lineare Regression eher ungeeignet.

# Plot with Plotly Express
px.scatter(x=train_dataset['bedrooms_per_rooms'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')

px.scatter(x=train_dataset['people_per_household'], y=train_dataset['median_house_value'], opacity=0.65, 
                trendline='ols', trendline_color_override='darkred')

sns.lmplot(data=train_dataset, x='people_per_household', y='median_house_value')

<seaborn.axisgrid.FacetGrid at 0x1a1836ebca0>

bedrooms_per_rooms und people_per_household zeigen einen negativen linearen Zusammenhang. Dies ists gut an der Linie im Plot verdeutlicht erkennbar. Bereits in der Korrelationsmatrix war der negative Zusammenhang dieser beiden Variablen zu median_house_value sichtbar. Die Variablen werden nachfolgend in die linearen Regressionsmodelle aufgenommen.

Modellierung¶

Statsmodel¶

Model 1 - OLS Regression¶

Als erstes Modell wird statsmodels OLS Regression verwendet.
Die sogenannten Filter Methode ist eine von drei Möglichkeiten bei der Feature Selektion. Gemäß dieser werden beispielsweise auf Basis einer Korrelationsmatrix die besten Feature ausgewählt und diese dem Modell zur Verfügung gestelltn.
Folgende Features wurden als geeignet identifiziert: median_income, ocean_proximity, bedrooms_per_rooms, people_per_household und eventuell housing_median_age.
Im ersten Modell werden die vier zuerst genannten Features verwendet.

# Fit Model
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household' , data=train_dataset).fit()
# Full summary
lm.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.619
Model:	OLS	Adj. R-squared:	0.619
Method:	Least Squares	F-statistic:	5312.
Date:	Tue, 18 Jan 2022	Prob (F-statistic):	0.00
Time:	01:03:32	Log-Likelihood:	-2.0580e+05
No. Observations:	16340	AIC:	4.116e+05
Df Residuals:	16334	BIC:	4.116e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	2.285e+04	4519.014	5.056	0.000	1.4e+04	3.17e+04
ocean_proximity[T.INLAND]	-6.729e+04	1415.385	-47.545	0.000	-7.01e+04	-6.45e+04
ocean_proximity[T.NEAR COAST]	1.676e+04	1428.672	11.729	0.000	1.4e+04	1.96e+04
median_income	4.338e+04	406.094	106.829	0.000	4.26e+04	4.42e+04
bedrooms_per_rooms	3.031e+05	1.3e+04	23.248	0.000	2.78e+05	3.29e+05
people_per_household	-1.049e+04	435.520	-24.097	0.000	-1.13e+04	-9640.930

Omnibus:	4763.697	Durbin-Watson:	1.971
Prob(Omnibus):	0.000	Jarque-Bera (JB):	22109.483
Skew:	1.349	Prob(JB):	0.00
Kurtosis:	8.019	Cond. No.	129.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

lm.params

Intercept                         22847.458850
ocean_proximity[T.INLAND]        -67294.434300
ocean_proximity[T.NEAR COAST]     16757.013965
median_income                     43382.483758
bedrooms_per_rooms               303138.303742
people_per_household             -10494.596699
dtype: float64

Die summary() Methode berechnet die wichtigsten Kennzahlen des Modells.
R-squared ist eine Kennzahl zur Beurteilung der Anpassungsgüte eines Modells und nimmt einen Wert zwischen 0% und 100% an. Als Bezugsbasis wird der Durchschnitt verwendet.
lm erreicht einen R-squared Wert in Höhe von 63,2%. Somit können 63,2% der median_house_value Streuung
In statistischen Methoden wird meist lieber der adjusted R-squared genutzt. Da dieser die degrees of freedom miteinbezieht, lässt sich durch adj. R-squared einen besseren Rückschluss auf die Gesamtpopulation der Datensätze ziehen. Dieser ist meistens etwas schlechter als R-squared. In diesem Modell sind beide Werte gleich hoch.
Mithilfe der Kennzahlen AIC und BIC soll ein overfitting des Modells verhindert werden. Für jedes zusätzliche Feature wird eine “Bestrafung” in die Kennzahlen miteinberechnet. Je kleiner die Werte, desto besser. In diesem Modell sind die Werte jedoch relativ hoch.
F-statistics ist die Menge der systematischen Varianz (MSM) geteilt durch die Menge der unsystematischen Varianz (MSR). Die Kennzahl gibt an, in welcher Höhe das Modell die Ergebnisausgabe der abhängigen Variable verbessert hat, verglichen zu Ungenauigkeit im Modell. Je höher F-statistics, desto besser. Ein F-statistics Wert in Höhe von 5618, wie in Modell lm, ist schon sehr gut.
Der Y-Achsenschnittpunkt des Modells liegt bei 61394,65. Die Steigung der einzelnen Features ist in der Spalte coef abzulesen. Zu interpretieren sind diese folgendermaßen: Vergrößert sich der Wert in Variable median_income um 1 und alle anderen Werte bleiben konstant, so ist das Ergebnis um 42189,93 größer als zuvor.
Die p-values der Features liegen alle bei 0,00. Die Nullhypothese (ß=0 => Feature hat keinen Einfluss auf die abhängige Variable) kann demnach für alle Features verworfen werden.
Die p-values für Jarque-Bera und Omnibus liegen ebenfalls bei 0,00. Die Nullhypothese, die Fehler seien normalverteilt, wird in diesem Falle verworfen. Es wird also angenommen, die Residuals seien nicht normalverteilt. Bei einer Datensetgröße > 50, kann dies aufgrund des Central Limit Theorems jedoch vernachlässigt werden.
Der Durbin-Watson Wert liegt bei knapp 2. Somit kann angenommen werden, dass die Fehler des Modells keine Korrelation aufweisen.

Regression Diagnostics¶

fig = sm.graphics.influence_plot(lm, criterion="cooks")
fig.tight_layout(pad=1.0)

Mithilfe des Influence Plots werden mögliche Ausreiser als blaue Punkte in der Größe entsprechend ihres Einflusses identifiziert. Dies können Werte mit hoher Abweichung in ihrem Y-Wert (outliers) oder X-Wert (high leverage) sein. In einer linearen Regression können solche Ausreise einen großen Einfluss auf das Gesamtmodell haben und somit die Performance des gesamten Modells erheblich verschlechtern. Als Faustregel werden Datenpunkte mit einer Cooks Distance in Höhe von 4/n (wobei n die Anzahl der Datensätze ist) als mögliche Ausreiser identifiziert.

# obtain Cook's distance 
lm_cooksd = lm.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])

n =  16340
Critical Cooks distance: 0.0002447980416156671
Int64Index([15778,  6776, 11511,  7263, 15723, 14756, 15687,  4638, 14806,
            15815,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=921) 
 [1.48759293e-03 2.66836227e-04 4.78884475e-04 5.45180812e-04
82637029e-04 5.55166261e-03 5.93629188e-04 9.64446395e-04
14112179e-04 8.23043391e-04 1.37179384e-03 9.31198962e-04
28353433e-04 3.14064566e-04 3.70703445e-04 7.92176597e-04
76686718e-04 2.53437155e-04 5.35087643e-04 4.34543876e-04
56743264e-04 4.37260421e-04 3.15946778e-04 3.85848620e-04
78654618e-04 9.85584673e-04 4.69834930e-02 5.19188525e-04
33272705e-04 2.46198971e-04 3.62984273e-04 3.72530337e-04
66246300e-04 8.19396201e-04 4.87217083e-04 1.63877589e-03
18039096e-03 1.10631997e-03 5.10220751e-04 2.98320853e-04
98424061e-04 2.83831187e-03 4.23210697e-04 4.84462529e-04
35832817e-04 2.56077726e-03 6.60015026e-04 3.00979233e-04
89758021e-04 9.03255145e-04 3.76539069e-04 4.30606581e-04
28854120e-04 3.86905704e-04 7.61360057e-04 3.48755331e-04
90570754e-04 5.32416641e-04 4.53912402e-04 5.91744313e-04
61448425e-04 4.27661004e-04 4.01913126e-04 7.92934770e-04
08537643e-04 4.60726908e-04 1.17463913e-03 2.51060698e-04
80930856e-04 5.12470023e-04 5.65809255e-04 3.27401509e-04
85242374e-04 1.21787762e-03 7.24325239e-04 2.92813642e-04
73003816e-04 4.68626337e-04 1.16110752e-03 4.63062325e-04
25016890e-04 4.83688775e-04 5.35658239e-04 1.94858180e-03
89813290e-04 8.85370793e-04 7.17970203e-04 6.42891705e-04
14955277e-04 1.02215153e-03 3.56207358e-04 8.67510310e-04
10491710e-03 7.99219333e-04 2.62696017e-04 2.88666981e-04
07854509e-04 8.83931891e-04 2.45081661e-04 1.06449798e-03
70938773e-04 3.09018081e-04 6.57379410e-04 1.06179504e-03
02894187e-04 9.26555582e-04 1.01293834e-03 3.20886325e-04
31556876e-03 3.21853302e-04 3.36049381e-04 2.50133964e-04
13847455e-04 2.52642306e-04 3.13725249e-04 7.12813234e-04
91857937e-04 5.77017000e-04 6.85531264e-04 2.60496718e-04
85522112e-04 2.61791680e-04 3.35125004e-04 7.60609675e-04
47063429e-04 8.90184001e-04 9.82453636e-04 3.85061491e-04
98611688e-04 4.21507557e-04 5.80036377e-04 2.97849695e-04
48910258e-04 4.04532613e-04 3.12791855e-02 7.55680452e-04
20510714e-04 2.91560940e-04 3.88296383e-04 2.61226212e-04
36742340e-04 7.57230958e-04 6.60575314e-04 5.89682358e-04
53674156e-04 5.67379890e-04 4.75221026e-04 3.19645157e-04
78539449e-04 3.63618083e-04 2.74057139e-04 5.54655710e-04
54335893e-04 2.75933612e-04 4.70547153e-04 6.15449463e-04
62179019e-04 3.00882646e-04 4.13638589e-04 2.65252379e-04
58798052e-04 3.59247990e-04 1.25803553e-03 4.35828212e-03
77124480e-04 3.30386011e-04 6.10246492e-04 3.95647548e-04
18201721e-04 3.56188388e-04 2.78201241e-04 2.45863813e-04
02007922e-04 2.50010455e-04 6.14141780e-04 2.92218505e-04
53697899e-04 2.46267670e-04 2.67098403e-04 7.51479956e-04
81424778e-04 3.27290114e-04 8.80851341e-04 1.20674591e-03
03575822e-03 5.18664791e-04 3.89595116e-04 2.68812640e-04
59722638e-04 1.20517420e-03 6.24208641e-04 3.72779920e-03
51415384e-04 3.72957326e-04 3.05334929e-04 5.52136850e-04
98068484e-03 3.29411537e-04 3.40411907e-04 2.56222465e-04
31335952e-04 3.11886294e-03 2.61795401e-04 3.49696722e-04
74889698e-04 4.24423862e-04 5.33480402e-03 3.76297264e-03
83074019e-03 8.29506869e-04 1.36747265e-03 4.79332644e-04
70922385e-04 3.90781122e-04 2.64913644e-04 5.65879731e-04
98086264e-04 2.94669543e-04 2.63610289e-04 3.68587927e-04
99298783e-04 4.35823599e-04 9.11267981e-04 9.39988920e+00
39861014e-04 2.69052780e-04 3.61433050e-04 6.04443335e-04
66277640e-04 7.66984878e-04 3.10471732e-04 4.34022341e-02
09218138e-04 3.68185782e-04 2.18766380e-03 3.56981125e-04
31902936e-04 1.39521400e-03 3.20738655e-04 4.98230233e-04
77448218e-04 3.32692672e-04 2.94224674e-02 9.53393326e-04
06615406e-04 5.79773700e-04 3.51514909e-04 3.61553795e-04
22957130e-03 4.63021859e-04 2.69086838e-04 4.12697044e-04
51813451e-04 9.70081562e-04 2.55405842e-04 1.29836196e-03
14407763e-03 4.27009366e-04 8.02132198e-04 2.68103122e-04
88662584e-04 5.24588063e-04 2.47282051e-04 3.09316188e-04
85607781e-04 3.62042318e-04 2.69988391e-04 2.60804951e-04
95435073e-04 2.64265161e-03 2.73619549e-04 3.09846370e-04
08953714e-03 3.19858062e-03 3.61945601e-04 4.48200055e-04
46775381e-04 2.45696522e-04 3.29958181e-04 4.37556828e-04
94720981e-04 8.15616897e-04 2.94586664e-04 2.46671711e-04
67944042e-04 5.20079922e-04 2.61685051e-04 3.00405843e-03
95447880e-04 2.69710091e-04 5.84435890e-04 7.01102896e-04
06679338e-03 5.43870147e-03 3.74467022e-04 2.94450188e-04
53604644e-04 3.86371760e-04 3.37832122e-04 4.08310873e-04
57857913e-04 2.53449617e-04 4.35238135e-04 9.32109329e-04
19634636e-04 3.70736472e-03 6.60444604e-04 2.13403879e-03
47318394e-04 3.93509961e-03 1.71125362e-03 3.14016243e-04
64171619e-04 3.54606811e-04 3.03053219e-04 4.38959252e-04
86291326e-03 4.18578111e-04 6.23666239e-04 2.74507153e-04
55770133e-04 2.39766662e-03 5.29707607e-04 1.47048722e-03
72050204e-04 5.91643897e-03 6.45921796e-04 1.71379171e-03
29776306e-04 4.08631064e-04 1.03297005e-03 5.77032046e-04
76044202e-04 6.17879439e-04 4.20098365e-04 1.41320439e-03
63198371e-04 2.39553721e-03 3.56444674e-04 3.76120928e-04
75049105e-04 4.20190411e-03 3.52751919e-04 8.64347341e-04
14599731e-03 5.21634440e-04 5.83742875e-04 2.69107084e-04
94255954e-04 6.08876932e-04 8.17456580e-04 4.13993717e-04
67213098e-04 7.19449438e-04 2.52664572e-04 3.51423970e-04
71139725e-04 4.30278077e-04 1.10036644e-03 4.48525821e-04
97875784e-04 2.53128372e-04 6.71543161e-04 2.57889557e-04
37631343e-04 5.31753783e-04 4.36362942e-04 3.48171741e-04
33669610e-03 6.27844689e-04 3.30353881e-04 2.63316932e-04
93430902e-03 3.84002969e-04 3.10408934e-04 8.02639776e-04
60331666e-04 3.55062676e-04 2.93513559e-04 7.77477806e-04
78901271e-04 2.56911171e-04 2.62926487e-04 2.66958932e-04
18433088e-04 4.64211041e-04 3.30739046e-04 2.52024154e-04
27278215e-03 2.84426835e-04 7.04396718e-04 6.33891174e-04
39793952e-04 2.99236648e-04 4.45190246e-04 5.35835774e-04
22579904e-04 1.14102901e-03 1.00360134e-03 3.28224562e-04
79656610e-04 3.27410721e-04 7.39804588e-04 3.89886247e-04
33262252e-04 3.94098964e-04 4.18819969e-04 8.19442679e-04
89986966e-04 5.36623487e-04 3.96495856e-04 2.53190918e-03
24780554e-04 2.76551422e-03 2.68771101e-04 3.54448329e-04
13845602e-03 2.46529930e-04 4.32895348e-03 4.15819743e-04
93225300e-04 3.82160961e-03 5.93387422e-04 6.77249628e-04
68111655e-04 1.48386533e-03 2.70829420e-04 2.60996883e-03
80995269e-04 5.33430607e-04 5.13911021e-04 2.56822250e-04
00464300e-04 5.89936555e-04 1.00222242e-03 4.49752787e-04
85344597e-03 4.42173403e-04 3.74656557e-04 4.05218668e-04
49437119e-04 8.48112319e-04 9.52066468e-04 3.17766109e-04
21040279e-04 2.67941768e-04 2.16147046e-03 5.23891516e-04
89009713e-03 5.82745086e-04 6.78943360e-04 3.07943494e-04
01012970e-04 4.48561295e-04 1.28186822e-03 8.61789035e-04
34727316e-04 2.73804253e-04 1.42202399e-01 3.50860752e-04
77507369e-03 4.14111841e-04 4.34040630e-04 8.47409581e-04
73707382e-04 9.36056215e-04 4.39638949e-04 2.47091843e-04
94815965e-03 3.41359419e-04 2.58914669e-04 1.29921590e-03
43443742e-04 7.08910616e-04 2.51070091e-04 7.51519537e-04
61757210e-04 2.28697927e-03 3.46303378e-04 2.72064007e-04
02842045e-04 1.08821274e-03 2.48982171e-03 2.53807767e-04
76891853e-04 2.47227261e-04 3.45705294e-04 4.99565882e-04
35401915e-03 2.70401709e-04 1.09188172e-02 2.91184063e-04
50205439e-04 3.35361471e-04 7.99127506e-04 4.69299374e-04
63469833e-04 7.40940278e-04 2.64031971e-04 3.29470784e-04
47023459e-03 3.40530222e-04 4.65887402e-03 4.50849207e-04
52524725e-04 8.97162356e-03 5.70524019e-04 7.51568069e-04
50790326e-04 2.85441357e-03 6.78405443e-03 3.43573733e-04
36485914e-04 2.98686005e-04 8.58965733e-04 2.54924272e-04
30030938e-04 2.56453063e-03 3.93964524e-04 3.95499423e-04
79750484e-04 5.50201808e-04 1.71339260e-03 3.37091835e-04
40192757e-04 8.24478158e-04 3.12070719e-04 1.35764062e-03
43487813e-04 3.12996792e-04 5.68496869e-04 3.98348038e-04
46630274e-04 4.95629626e-04 2.57410000e-04 8.08578807e-04
64273317e-04 2.82160955e-04 3.65064309e-04 2.70263308e-04
45828059e-04 8.38319883e-04 3.25170223e-04 3.06324121e-04
33351819e-04 6.77037097e-03 2.75609130e-04 1.01649062e-03
04145849e-04 5.15573187e-04 3.88917919e-04 1.46058927e-03
71213447e-04 1.08778812e-03 2.56939691e-04 6.19395605e-04
94834315e-04 4.06971618e-04 2.78217893e-03 3.23101550e-04
74127917e-04 1.76364447e-03 8.36296107e-04 3.80469230e-04
70007674e-04 1.17371929e-03 6.95576355e-04 3.32710828e-04
59662126e-04 5.69312329e-04 9.26900479e-04 8.42265531e-04
10165863e-04 2.70328570e-03 7.71394059e-04 2.68445474e-04
11973612e-04 3.32461369e-04 9.80372135e-04 8.81917284e-04
55373959e-04 3.29644215e-04 6.21269904e-04 3.76872079e-04
50203079e-04 4.52833892e-04 2.97099958e-04 5.22922148e-04
92262851e-04 2.20052284e-03 3.83532304e-04 6.32149739e-04
65631266e-04 4.72368031e-04 7.89694271e-04 8.37232630e-04
74232894e-04 5.64497691e-04 2.94122337e-03 2.63917959e-04
04321725e-04 2.89069692e-04 6.90659396e-04 2.63585485e-04
91107222e-04 2.83070273e-04 9.32592868e-04 5.66913333e-04
08362849e-04 3.36558018e-03 2.76066020e-03 2.87497857e-04
19049137e-03 4.52875210e-04 2.67604795e-04 4.02825302e-04
05156479e-04 3.06836938e-04 2.78851619e-04 1.74883211e-03
64067879e-04 6.07388648e-04 7.58277615e-04 7.57273054e-04
29332460e-04 2.81986374e-04 2.82743794e-04 4.05569007e-04
83592900e-04 7.04332174e-04 2.75563066e-04 4.31355869e-04
42656246e-04 1.82385222e-03 5.22656398e-04 3.51338107e-04
79176372e-04 1.03947054e-03 1.07116863e-03 1.48117322e-03
65401165e-04 4.37486711e-04 1.09231708e-03 2.45871300e-04
31703787e-04 3.31251709e-04 4.01178501e-04 1.48896621e-03
56162782e+00 5.61537488e-04 4.83215948e-04 1.21070439e-03
75200724e-04 4.01127609e-04 6.27011828e-04 2.66717361e-04
95563009e-04 3.02506345e-04 1.44876067e-03 3.00789351e-04
91476001e-04 4.39010114e-04 3.17796695e-04 2.98886372e-04
91550120e-04 4.73051682e-04 6.31527738e-04 9.31808401e-04
03824501e-04 3.79034197e-04 2.59611879e-04 6.49811471e-04
87937183e-04 3.80029796e-03 5.59943324e-04 3.17157107e-04
92633169e-04 4.14709560e-04 6.86773617e-01 2.48380968e-03
90093477e-04 4.36361431e-04 2.59591922e-04 2.47050144e-04
59685093e-04 3.01684047e-04 3.02507235e-04 5.18916384e-04
84954521e-04 4.74586441e-04 3.04089487e-04 4.26933571e-04
91491532e-03 3.17801973e-04 2.47929448e-04 9.06269529e-04
94183033e-03 3.48312422e-04 5.50505354e-04 4.58721438e-04
01822960e-03 2.52072586e-04 4.91210068e-04 5.05974483e-04
22862009e-03 4.23452403e-04 1.21258864e-03 4.38637903e-04
89214946e-04 3.75037377e-04 2.58273666e-04 2.52480169e-04
55347542e-04 5.32865237e-04 1.62222644e-03 2.96881492e-04
81489423e-04 5.10193227e-04 2.53749020e-04 2.53913870e-04
48709098e-04 2.63204276e-04 3.58335313e-04 5.10050913e-04
06070501e-04 2.69330637e-04 2.98364767e-04 3.28267579e-04
54259807e-04 2.46176219e-04 3.72063464e-04 3.41298513e-04
07934920e-04 3.66198609e-04 3.66939372e-03 3.67570572e-03
90089084e-04 3.55853552e-04 2.60724438e-04 4.37018363e-04
57131144e-04 5.06057744e-04 2.78888173e-04 4.97652487e-04
39670036e-04 8.05463091e-04 3.51184750e-04 2.62145544e-04
74375146e-04 3.82508222e-04 3.68841183e-03 4.19717652e-04
55690629e-04 2.99726571e-04 2.89178308e-04 5.08179827e-04
03202235e-04 2.94196485e-04 2.81584302e-04 9.80940709e-03
16255126e-04 2.47946510e-04 2.63466143e-03 2.88074200e-04
67563020e-04 3.03789671e-04 5.36379052e-04 4.62299803e-04
80226710e-04 3.77007575e-04 7.47209568e-04 2.66495405e-04
13772364e-04 7.21450302e-04 5.47056613e-04 7.37967987e-04
27841466e-04 3.76690865e-04 4.37975867e-04 3.09631795e-03
84687078e-04 2.55618180e-04 3.45388186e-03 2.93628277e-04
61960005e-04 3.94595188e-04 5.89990377e-03 9.02572938e-04
82159330e-04 4.91128582e-04 7.25198188e-04 3.00210510e-04
85064819e-04 6.80836379e-04 6.44973092e-04 3.17077777e-04
22945658e-04 4.01621666e-04 1.20344371e-03 1.07300818e-02
66737814e-04 3.12277612e-04 3.20768810e-04 2.68265320e-04
91625279e-03 2.79387537e-04 4.52872521e-04 3.01871767e-04
26111818e-04 4.64855935e-04 3.62488941e-04 3.18872152e-04
99944069e-04 4.51085526e-04 4.98276715e-04 3.30513367e-04
31357117e-04 2.45649148e-04 3.97569437e-04 2.53923270e-04
79669582e-04 2.97081472e-04 2.65974062e-04 2.83734175e-04
45992907e-04 6.82390231e-04 3.25452929e-04 4.83621995e-04
95153421e-04 2.82123615e-03 2.37722866e-03 2.60311484e-04
78938163e-04 2.00983379e-03 9.78866811e-04 2.89492097e-04
91580342e-04 3.87343121e-04 6.97254658e-04 4.36301857e-04
17114784e-03 3.22186946e-04 3.85647265e-04 2.79761337e-03
39186665e-04 4.94425448e-04 2.62787157e-04 5.30634894e-04
94923616e-04 2.53532330e-04 7.41679657e-04 5.64260677e-04
77355811e-04 7.28817918e-04 3.09293347e-04 2.98152988e-04
19426368e-04 2.69393681e-04 2.93032510e-04 6.93512946e-04
26959929e-04 4.69588837e-04 1.08769402e-03 5.98038443e-04
76007963e-04 2.87637693e-04 4.54607469e-04 3.53033837e-04
58269976e-03 3.65366914e-04 3.06447575e-03 2.58277589e-04
63480494e-04 2.81210750e-04 6.35883420e-04 3.11540223e-04
52162584e-04 3.06498378e-04 3.74048116e-04 2.54321920e-04
80908175e-04 4.24951356e-03 2.84954354e-04 4.18691693e-01
29198280e-04 2.34973565e-03 4.18652522e-04 4.86206706e-04
26277968e-03 2.60487497e-04 1.23711073e-03 5.57395720e-04
58085248e-04 3.92749045e-04 3.82798861e-04 2.85589944e-04
96116397e-03 7.62436805e-04 3.62626773e-04 3.44899705e-04
37202800e-04 2.92916378e-04 2.49218239e-04 7.95948605e-04
12065059e-04 4.58312754e-04 4.16739905e-04 5.51838715e-04
87625848e-04 4.50237635e-04 2.06161193e-03 3.32015810e-04
83970721e-04 2.89200280e-04 3.67273098e-04 2.84789133e-04
30385790e-04 2.72907722e-04 3.10285352e-04 9.78476629e-04
03511525e-03]

In diesem Modell werden 963 Datensätze mit einer Cooks Distance > 4/n als mögliche Ausreiser identifiziert. Beispielsweise könnte versucht werden, den Einfluss dieser Datensätze auf das Modell durch Datentransformation (bspw. durch Verwendung des Logarithmus) zu verringern. Allerdings hilft diese Methode der Modellperformance genau so oft, wie sie ihr schadet. Daher wird hier auf eine Datentransformation verzichtet und die Datenpunkte stattdessen aus dem Trainings Set entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 15419 entries, 14185 to 17879
Data columns (total 14 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   longitude             15419 non-null  float64 
 1   latitude              15419 non-null  float64 
 2   housing_median_age    15419 non-null  int64   
 3   total_rooms           15419 non-null  int64   
 4   total_bedrooms        15419 non-null  int64   
 5   population            15419 non-null  int64   
 6   households            15419 non-null  int64   
 7   median_income         15419 non-null  float64 
 8   median_house_value    15419 non-null  float64 
 9   ocean_proximity       15419 non-null  category
 10  people_per_household  15419 non-null  float64 
 11  rooms_per_household   15419 non-null  float64 
 12  bedrooms_per_people   15419 non-null  float64 
 13  bedrooms_per_rooms    15419 non-null  float64 
dtypes: category(1), float64(8), int64(5)
memory usage: 1.7 MB

fig = sns.scatterplot(data=train_dataset2, x='median_income', y='median_house_value')

Im Vergleich zum ursprünglichen Trainings Set, welches im Kapitel EDA ausführlich analysiert wurde, wurden einige X- und Y-Ausreiser entfernt. Dadurch hat sich die waagrechte Punktelinie bei 500000$ medin_house_value verkleinert und Ausreiser mit hohen X- und niedrigen Y-Werten wurden entfernt.

fig = sns.displot(data=train_dataset2, x='median_house_value')

Durch die Entfernung der Ausreiser hat sich das Histogramm von median_house_value etwas weiter normalisiert.

# Fit Model without outliers
lm = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household', data=train_dataset2).fit()
lm.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.739
Model:	OLS	Adj. R-squared:	0.739
Method:	Least Squares	F-statistic:	8725.
Date:	Tue, 18 Jan 2022	Prob (F-statistic):	0.00
Time:	01:04:55	Log-Likelihood:	-1.8952e+05
No. Observations:	15419	AIC:	3.790e+05
Df Residuals:	15413	BIC:	3.791e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	3.341e+04	4404.814	7.585	0.000	2.48e+04	4.2e+04
ocean_proximity[T.INLAND]	-6.1e+04	1098.946	-55.509	0.000	-6.32e+04	-5.88e+04
ocean_proximity[T.NEAR COAST]	6838.7281	1120.708	6.102	0.000	4642.008	9035.448
median_income	4.689e+04	365.489	128.285	0.000	4.62e+04	4.76e+04
bedrooms_per_rooms	3.58e+05	1.21e+04	29.558	0.000	3.34e+05	3.82e+05
people_per_household	-2.492e+04	583.778	-42.686	0.000	-2.61e+04	-2.38e+04

Omnibus:	958.072	Durbin-Watson:	1.986
Prob(Omnibus):	0.000	Jarque-Bera (JB):	1206.262
Skew:	0.599	Prob(JB):	1.16e-262
Kurtosis:	3.664	Cond. No.	154.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Durch die Entfernung der Ausreiser haben sich R-squared, adjusted R-squared und F-statistics deutlich verbessert.
AIC und BIC sind auch etwas kleiner geworden und haben sich somit verbessert.
Die p-values der Features bleiben unverändert.
Auch an den p-values für Jarque-Bera und Omnibus hat sich kaum etwas verändert.
Durbin-Watson liegt nun bei einem Wert <> 2. Generell ist bei einem Wert von > 2 eine negative Korrelation der Modellfehler anzunehmen.

fig = sm.graphics.plot_partregress_grid(lm)
fig.tight_layout(pad=0.2)

Die Partial Regression Plots visualisieren die Relation einer unabhängigen Variablen in Verbindung mit den anderen abhängigen Variablen im Modell in Bezug auf die abhängige Variable.
Die Plots zeigen eine negative Abhängigkeit von der Ausprägung I.INLAND von ocean_proximity und peope_per_household. Median_income und bedrooms_per_rooms scheinen eine positive Abhängigkeit auf median_house_value zu haben. Gemäß dem Diagramm hat die Variable ocean_proximity[T.NEAR COAST] keinen Einfluss auf median_income. Dies steht im Gegensatz zum p-value 0,00 dieser Variablen, nach dessen Interpretation die Variable einen Einfluss auf die abhängige Variable hat.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)

Auch diese Plots lassen darauf schließen, dass die Variable _ocean_proximity[T.NEAR COAST] das Modell nicht verbessert. Sowohl Partial Regression Plot als auch CCPR Plot visualisieren eine Nicht-Linearität zwischen der Variablen und median_house_value.
Die Residuals scheinen dagegen gleichmäßig verteilt zu sein.

fig = sm.graphics.plot_regress_exog(lm, "ocean_proximity[T.INLAND]")
fig.tight_layout(pad=1.0)

Die negative Korrelation ist sowohl im Partial Regression Plot als auch im CCPR Plot zu erkennen. Die Residuals sind gleichmäßig verteilt und lassen daher auf Homoskedastizität rückschließen. Das Diagramm zu Y and Fitted vs. X ist aufgrund der vielen Datenpunkte schwer zu analysieren.

fig = sm.graphics.plot_regress_exog(lm, "people_per_household")
fig.tight_layout(pad=1.0)

Auch hier ist eine negative Abhängigkeit der abhängigen zur unabhängigen Variablen zu erkennen. Linearität scheint also vorzuliegen. Allerdings ist im Residual Plot keine konstante Varianz der Fehlerstreuung zu erkennen. Mit steigendem X-Wert verringert sich die Fehlerstreuung und liegt nur noch im positiven Bereich. Mit zunehmender people_per_household wird tendenziell ein höherer median_house_value mit dem Modell hervorgesagt, als tatsächlich vorliegt. Dies spiegelt sich auch im Plot zu Y and Fitted vs. X wieder. Bei dieser Variablen liegt somit Heteroskedastizität vor. Die Ursache hierfür könnte darin liegen, dass im höheren Wertebereich von people_per_household deutlich weniger Datensätze vorliegen als im niedrigeren Bereich und das Modell somit nicht ausreichend für diese Wertebereiche trainiert wurde. Womöglich bildet auch eine polynomiale Modellierung für diese Variable die Relation der Variablen besser ab.

fig = sm.graphics.plot_regress_exog(lm, "median_income")
fig.tight_layout(pad=1.0)

Für median_income ist die Linearität im Partial Regression Plot und CCPR Plot deutlich zu erkennen. Jedoch nimmt auch hier die Streuung der Residuals im höheren X-Wertebereich ab.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR COAST': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST' }

train_dataset.ocean_proximity = train_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

Model 2 - OLS Regression¶

Nachfolgend wird ein weiteres Modell mit zusätzlicher Verwendung der Variablen housing_median_age trainiert.

# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset).fit()
lm2.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.627
Model:	OLS	Adj. R-squared:	0.627
Method:	Least Squares	F-statistic:	5486.
Date:	Tue, 18 Jan 2022	Prob (F-statistic):	0.00
Time:	01:08:50	Log-Likelihood:	-2.0563e+05
No. Observations:	16340	AIC:	4.113e+05
Df Residuals:	16334	BIC:	4.113e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	-7.016e+04	4213.120	-16.652	0.000	-7.84e+04	-6.19e+04
ocean_proximity[T.NEAR COAST]	6.563e+04	1348.977	48.655	0.000	6.3e+04	6.83e+04
median_income	4.445e+04	406.063	109.460	0.000	4.37e+04	4.52e+04
bedrooms_per_rooms	2.998e+05	1.29e+04	23.243	0.000	2.74e+05	3.25e+05
people_per_household	-1.096e+04	428.532	-25.570	0.000	-1.18e+04	-1.01e+04
housing_median_age	1005.0830	46.230	21.741	0.000	914.467	1095.699

Omnibus:	4984.528	Durbin-Watson:	1.968
Prob(Omnibus):	0.000	Jarque-Bera (JB):	27788.201
Skew:	1.356	Prob(JB):	0.00
Kurtosis:	8.785	Cond. No.	766.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

R-squared und adj. R-squared haben sich im Vergleich zum ersten Modell verbessert.
Auch F-statistics hat sich durch die Entfernung der Ausprägung von ocean_proximity und Aufnahme des Features median_housing_age verbessert.
Alle übrigen Kennzahlen haben sich kaum verändert.

# obtain Cook's distance 
lm_cooksd = lm2.get_influence().cooks_distance[0]

# get length of df to obtain n
n = len(train_dataset["median_income"])
print('n = ', n)

# calculate critical d
critical_d = 4/n
print('Critical Cooks distance:', critical_d)

# identification of potential outliers with leverage
out_d = lm_cooksd > critical_d

# output potential outliers with leverage
print(train_dataset.index[out_d], "\n", 
    lm_cooksd[out_d])

n =  16340
Critical Cooks distance: 0.0002447980416156671
Int64Index([15778,  6776, 10005, 11511,  7263, 15723, 14756, 15687,  4638,
            14806,
            ...
            15809,  8866,   199, 18498,    59, 10749,  8857, 13154,  5470,
            16642],
           dtype='int64', length=834) 
 [1.34620956e-03 3.28749861e-04 2.57215207e-04 3.35039575e-04
51518260e-04 3.61919469e-04 7.91550903e-03 6.11015995e-04
49112668e-04 4.17870132e-04 8.15554753e-04 1.30731202e-03
20416904e-03 3.73036827e-04 5.62180943e-04 7.70163204e-04
46624084e-04 2.53753610e-04 3.01304709e-04 5.38305099e-04
35811574e-04 2.54781398e-04 4.29838436e-04 3.35159404e-04
69843859e-04 3.46408581e-04 2.57119718e-04 5.47776435e-02
46822263e-04 4.27225287e-04 3.34569001e-04 8.73800145e-04
12264946e-04 1.65774076e-03 2.06412271e-03 1.04840879e-03
72011719e-04 5.81620070e-04 5.20076218e-04 5.13152695e-03
17172666e-04 9.28482089e-04 4.71110244e-03 3.95928870e-04
95769203e-04 2.62218028e-04 9.37446313e-04 2.83036311e-04
06475831e-04 4.26676622e-04 3.01809213e-04 5.76534941e-04
38711782e-04 3.70526120e-04 4.21168257e-04 1.19851082e-03
09081908e-04 6.12841623e-04 3.59409818e-04 3.20014681e-04
94920532e-04 9.53529426e-04 2.57227742e-04 8.06320199e-04
55183702e-04 4.39410335e-04 6.62378271e-04 5.62770048e-04
31933928e-04 7.31914431e-04 1.00363583e-03 5.29889822e-04
33576686e-04 4.95675084e-04 8.89908119e-04 2.83041383e-04
55849228e-04 2.16001576e-03 3.09497542e-04 6.81234101e-04
08301229e-04 6.74865681e-04 3.14174140e-04 9.77863893e-04
98076323e-04 8.43691811e-04 1.18174422e-03 9.89622636e-04
13546336e-04 9.31650342e-04 1.20073880e-03 3.56089239e-04
11493487e-04 4.13533793e-04 1.64662507e-03 3.81249929e-04
84393006e-04 6.50545586e-04 2.18054214e-03 2.78936483e-04
52932373e-04 2.52777940e-04 2.66069166e-04 5.83492758e-04
39562190e-04 6.77350435e-04 6.48707359e-04 4.21307444e-04
91646476e-04 2.47710633e-04 2.53301281e-04 3.56535739e-04
43918045e-04 3.77300968e-04 8.44956605e-04 1.30709420e-03
85992118e-04 4.39369242e-04 3.80927914e-04 3.82280250e-04
09153132e-04 3.44268183e-04 5.46423976e-04 4.20173053e-04
28627490e-02 7.35095993e-04 4.24027988e-04 7.22312913e-04
32786505e-03 6.99537001e-04 2.50626133e-04 5.46342106e-04
94243518e-04 3.65875470e-04 3.13187264e-04 2.70271254e-04
25703718e-04 6.93446296e-04 2.65491825e-04 3.49024635e-04
08335882e-04 4.05092100e-04 2.45889472e-04 4.39110074e-04
85502846e-04 9.38717658e-04 1.34362145e-03 4.80589164e-03
98106969e-04 3.22240865e-04 3.15091984e-04 2.56911300e-04
46225434e-04 4.22497484e-04 3.46518228e-04 6.70520360e-04
22863087e-03 7.43148088e-04 2.54542565e-04 2.56665208e-04
40384851e-03 3.28558370e-03 3.11047780e-04 3.43185408e-04
43474503e-03 7.15000996e-04 4.19995923e-03 7.80899608e-04
75125809e-04 2.79783608e-04 3.39099740e-04 1.19937066e-03
88403117e-03 3.65830328e-04 5.86482204e-04 2.56515079e-04
23646462e-03 2.54733684e-04 4.59159574e-04 8.87961889e-03
74045810e-03 4.00718231e-03 5.46171186e-04 1.50078936e-03
56173603e-04 4.34240443e-04 2.61806316e-04 3.31107538e-04
18033323e-03 2.47950011e-04 3.32658507e-04 3.22346407e-04
61654415e-04 6.16450930e-04 9.91135344e-04 1.06921629e+01
82232779e-04 3.13847258e-04 3.09257186e-04 3.06408629e-04
56403994e-04 1.05400214e-03 4.00527883e-02 2.84546161e-04
53531661e-04 2.11139733e-03 3.83398773e-04 4.32087727e-04
45600603e-03 6.41920570e-04 6.14620185e-04 3.14498312e-04
97336263e-04 3.07533280e-02 8.85053202e-04 4.54537558e-04
20079773e-04 3.73657618e-04 2.85897551e-04 1.45010004e-03
47302087e-04 2.86417016e-04 3.86170479e-04 3.64811817e-04
04250746e-04 8.96897085e-04 3.31137357e-03 4.58857362e-04
13777854e-04 2.81677164e-04 3.09589735e-04 5.65060691e-04
58478200e-04 3.69278949e-04 1.84718845e-03 2.58189555e-04
25494577e-04 4.15475232e-04 4.83547910e-03 2.45981050e-04
75498461e-04 7.50737541e-04 4.33679600e-03 3.02966513e-04
16482267e-03 4.43903600e-04 4.39257693e-04 5.23686883e-04
43692280e-04 6.87454121e-04 3.78535731e-04 2.49192273e-04
34747272e-04 3.39985468e-04 3.30260839e-04 2.66248512e-03
09325047e-04 3.81360279e-04 6.38501803e-04 7.93860341e-04
17777251e-03 4.86447450e-03 4.16377804e-04 3.53579296e-04
94406173e-04 2.61202267e-04 3.66836210e-04 2.62061270e-04
38031251e-04 4.27777058e-04 9.13378989e-04 3.41456845e-04
64618565e-03 4.72926716e-04 1.77649059e-03 5.02697205e-04
67504825e-03 1.71681800e-03 5.55103694e-04 2.91051755e-04
55591228e-04 4.59132619e-04 2.07003198e-03 3.97897537e-04
05201888e-04 6.45379090e-04 3.27399025e-04 4.93574804e-04
46207746e-03 1.06687385e-03 1.27664267e-03 5.98190815e-04
69362545e-03 6.47563661e-04 1.41449488e-03 3.35907072e-04
31922192e-04 9.62879839e-04 5.99119839e-04 5.94298294e-04
72450549e-04 1.33284800e-03 2.22449093e-03 2.94392265e-04
60548234e-04 2.46831910e-04 5.21180548e-03 3.74112820e-04
00529729e-04 1.83317278e-03 2.10050482e-03 2.46498466e-04
52625587e-04 5.69750888e-04 2.70235120e-04 4.73789257e-04
37396918e-04 1.33110610e-03 3.97863971e-04 3.91811176e-04
56690305e-04 3.49347749e-04 3.60971784e-04 2.55520221e-04
55568409e-04 3.07803326e-04 1.03831598e-03 2.82099059e-04
35719194e-04 2.63385336e-04 7.30688902e-04 2.46434487e-04
32206636e-04 5.71495033e-04 4.55220555e-04 4.58525682e-04
51076913e-04 7.81748699e-03 6.88520513e-04 4.47783205e-04
38917952e-03 4.58028391e-04 8.53315502e-04 7.98043674e-04
06666704e-04 3.57944205e-04 2.97907283e-04 6.60775311e-04
15356513e-04 2.64416584e-04 4.01810556e-04 4.11109209e-04
17853473e-03 3.03560371e-04 7.12310609e-04 7.81473143e-04
45443831e-04 4.31752702e-04 3.55261927e-04 4.28943588e-04
18408870e-03 1.26262767e-03 3.30624525e-04 3.11377904e-04
37247025e-04 4.34417931e-04 3.13091517e-04 3.86335627e-04
27376743e-04 1.69310467e-03 3.75865197e-04 2.94169921e-04
09415944e-03 5.42417148e-04 5.02146968e-03 2.72579606e-04
40573948e-04 1.36489576e-03 3.93579498e-03 5.59166602e-04
11084962e-04 3.49829003e-04 3.02065944e-03 4.74636086e-04
17112634e-04 5.40381997e-04 2.21043242e-03 2.99376796e-04
78605704e-03 2.65136147e-04 3.62158887e-04 5.37258826e-04
16663410e-04 3.34015704e-04 2.37806674e-03 3.87874124e-04
60018714e-03 4.77441784e-04 3.16011510e-04 8.63052283e-04
69780131e-04 2.87389645e-04 4.19416247e-04 5.04048059e-04
01334240e-03 3.63296782e-04 3.86119046e-03 6.21399855e-04
91604511e-04 3.00672616e-04 3.07340296e-04 8.51472390e-04
39838656e-03 6.99686628e-04 3.88329899e-04 1.44844319e-01
06391695e-03 3.81700667e-04 4.36672272e-04 6.56546831e-04
04080380e-04 1.40540348e-03 4.24763158e-04 5.29649070e-03
49851273e-04 1.37361535e-03 5.45285648e-04 7.92267019e-04
89794507e-04 2.05090829e-03 3.00435326e-04 3.57210150e-04
07742500e-04 1.18352972e-03 2.43728565e-03 2.82272469e-04
89905770e-04 2.25566831e-03 3.03114176e-04 1.30068085e-02
58994479e-04 5.53658856e-04 6.12830952e-04 5.01723767e-04
12629182e-04 3.77644531e-04 1.55498141e-03 3.05627947e-04
94422552e-03 3.57403852e-04 8.78405940e-03 3.13205038e-04
27447728e-04 7.44573473e-04 5.15549175e-03 7.62643519e-03
32130683e-04 2.61203844e-04 1.29768456e-03 2.48558136e-04
25434452e-04 3.67663029e-03 4.20029578e-04 4.90839693e-04
73133315e-04 3.38691084e-04 1.89005180e-03 4.28670600e-04
55054527e-04 3.18588498e-04 7.39903481e-04 2.31834749e-03
67185876e-04 3.65229676e-04 5.75161805e-04 3.15816880e-04
74515715e-04 5.13067456e-04 6.89572009e-04 2.77145250e-04
70737869e-04 2.64762172e-04 5.27422472e-04 9.84707517e-04
43212824e-03 1.70314323e-03 4.62703809e-04 3.92885812e-04
20675602e-04 1.95409614e-03 2.88329437e-04 1.32460983e-03
43856512e-04 4.88842954e-04 2.69788415e-04 3.59022621e-03
48583713e-04 1.65698848e-03 2.23761619e-03 3.61044126e-04
20461757e-04 3.31346884e-04 2.44398609e-03 7.30058445e-04
79412548e-04 6.82860694e-04 8.29685275e-04 1.15423848e-03
92749800e-03 8.76630693e-04 2.63359256e-04 5.61268507e-04
13278888e-04 1.70939726e-03 2.71500955e-04 8.65533017e-04
66559196e-04 6.66489288e-04 3.82956960e-04 3.70377201e-04
62988149e-04 2.48974052e-04 2.90398752e-04 3.77658563e-04
49335881e-03 2.53154039e-04 2.66962782e-04 9.60055604e-04
36687000e-04 6.75338838e-04 2.07363805e-03 7.50365981e-04
87624234e-04 4.74190135e-03 2.87536148e-04 3.54576216e-04
61051300e-04 3.07458358e-04 3.37358573e-04 3.19957211e-04
92073680e-04 6.06543953e-04 5.92986873e-04 3.08890432e-04
43232872e-03 3.15607985e-03 3.87409151e-04 1.20192208e-03
33069594e-04 7.01263187e-04 3.66029598e-04 1.51065092e-03
10965707e-04 9.46617435e-04 1.33914359e-03 5.42132639e-04
55992635e-04 2.79273013e-04 4.33360234e-04 9.02970286e-04
22861701e-04 5.34084933e-04 4.66705558e-04 3.30393525e-03
95158501e-04 5.58204627e-04 2.76748391e-04 9.28388580e-04
60918152e-04 8.80922972e-04 1.07105882e-03 3.69483391e-04
64798647e-04 2.73824841e-04 9.66985466e-04 2.83190149e-04
82375269e-04 2.60887897e-04 3.25360840e-04 1.15909394e-03
02187581e+00 5.43475220e-04 4.02644115e-04 8.64339222e-04
18907455e-04 3.24992180e-04 7.38184531e-04 2.51052121e-04
42272514e-04 3.04544105e-04 3.63665615e-04 1.02912398e-03
87378521e-04 7.45509328e-04 3.26415591e-04 2.52550349e-04
40003285e-04 4.24117891e-04 4.90316626e-04 4.86412402e-04
40754977e-04 3.39882911e-04 6.81558602e-04 3.32647839e-04
38165023e-03 4.58021592e-04 6.59793292e-04 5.37127516e-04
46891438e-01 2.20109675e-03 3.91297892e-04 4.89698239e-04
60150085e-04 3.59658963e-04 3.11584443e-04 3.02795328e-04
18784806e-04 2.78422616e-04 4.22763551e-04 2.65366015e-04
65452108e-04 2.68691754e-04 2.98144020e-04 3.86472350e-03
38594356e-04 8.99398215e-04 5.10549736e-03 4.16733375e-04
03199774e-04 7.18913911e-04 4.79645149e-04 5.40162714e-03
87982015e-04 6.20419951e-04 4.25703669e-03 4.71147289e-04
04261716e-03 3.70599935e-04 4.27817891e-04 6.19133596e-04
53943173e-04 3.41160201e-04 1.66679303e-03 6.29287626e-04
17842252e-04 3.12862593e-04 7.27747807e-04 2.61667141e-04
64546973e-04 4.43050694e-04 2.89239694e-04 3.40418715e-04
36218958e-04 2.98524754e-04 8.09693242e-04 1.15337376e-03
63411919e-04 3.58669103e-03 2.98727727e-03 2.53558984e-04
42658703e-04 5.46202372e-04 3.09921088e-04 4.24199927e-04
66680028e-04 2.52166862e-04 4.51032150e-04 3.98717184e-04
78799300e-04 4.43084627e-04 2.66797596e-04 2.70016395e-04
52045620e-04 6.71406313e-04 3.51355939e-03 4.71552792e-04
67398242e-04 3.26641206e-04 4.89193724e-04 5.54156290e-04
55441606e-04 8.70449188e-03 1.01979909e-03 4.47454853e-03
67577330e-04 6.03933837e-04 4.19977124e-04 3.30883814e-04
59675243e-04 7.82055746e-04 5.70203443e-04 2.82550451e-04
75940421e-04 5.66976263e-04 7.23078868e-04 5.58070994e-04
80964154e-04 4.70836891e-04 4.63370304e-03 2.93061722e-04
82097527e-04 5.25394787e-04 4.39058608e-03 3.02814301e-04
93885720e-04 2.72975994e-04 6.30941339e-04 6.34256527e-03
93264396e-03 2.63895051e-04 3.90119304e-04 6.49958891e-04
39334191e-04 2.80016960e-04 6.70981971e-04 7.96887450e-04
61031562e-04 2.99531345e-04 3.92412896e-04 9.15393354e-04
65937206e-03 2.59496858e-04 3.04018109e-04 3.35152782e-04
47545315e-03 3.28849297e-04 3.09466337e-04 2.57128750e-04
79134179e-04 2.84444223e-04 7.04680802e-04 3.08216566e-04
19066173e-04 3.31656683e-04 4.99721752e-04 3.88382147e-04
62130992e-04 3.38338725e-04 3.33453243e-04 3.36783921e-04
53353487e-04 4.69560563e-04 1.01286041e-03 2.65157668e-04
44166808e-04 5.14610488e-04 3.25086702e-04 4.05530866e-03
61966032e-03 3.81661774e-04 1.26797417e-03 1.00604063e-03
47097398e-04 2.54264288e-04 3.76542325e-04 5.62680829e-04
39275303e-04 3.78306069e-03 4.09577971e-04 4.23492354e-03
67446416e-04 4.53267869e-04 2.61729406e-04 4.40048926e-04
43015952e-04 5.68428522e-04 2.68099048e-04 3.75476250e-04
99763254e-04 3.25127465e-04 2.44953335e-04 5.51778636e-04
10795171e-04 3.74624073e-04 6.56877859e-04 4.18323970e-04
79464330e-04 7.98927851e-04 3.65242322e-04 2.69035991e-04
39453306e-04 4.82734528e-04 3.59781628e-04 6.59773820e-03
12038525e-04 2.57552160e-04 3.11064867e-03 5.09224651e-04
60255237e-04 2.55022118e-04 5.33067322e-04 3.11642549e-04
90908224e-04 6.34810463e-04 4.38153846e-04 5.71802541e-03
88659684e-04 3.99323486e-01 5.10145404e-04 2.98228342e-03
12344844e-04 5.30124261e-04 9.11032421e-04 1.47777500e-03
02556616e-03 4.81595818e-04 4.69679318e-04 2.80467578e-04
75736314e-03 1.01866312e-03 4.80632918e-04 3.61741717e-04
65964061e-04 2.66946243e-04 1.07338565e-03 2.64243381e-04
62524321e-04 2.64787790e-04 3.70798360e-04 2.91009002e-04
70681350e-04 3.49234458e-04 3.54200151e-04 1.38738166e-03
21375424e-04 3.65064441e-04 3.17517938e-04 4.03136086e-04
73400173e-04 8.99485876e-04 2.73811960e-04 3.31912000e-04
43051451e-04 5.70143496e-03]

Auch in Modell lm2 werden mögliche outliers und Punkte mit high leverage identifiziert und aus dem Datenset entfernt.

#drop outliers identified by Cook's distance
train_dataset2 = train_dataset.drop(train_dataset.index[out_d])
train_dataset2.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 15506 entries, 14185 to 17879
Data columns (total 14 columns):
 #   Column                Non-Null Count  Dtype   
---  ------                --------------  -----   
 0   longitude             15506 non-null  float64 
 1   latitude              15506 non-null  float64 
 2   housing_median_age    15506 non-null  int64   
 3   total_rooms           15506 non-null  int64   
 4   total_bedrooms        15506 non-null  int64   
 5   population            15506 non-null  int64   
 6   households            15506 non-null  int64   
 7   median_income         15506 non-null  float64 
 8   median_house_value    15506 non-null  float64 
 9   ocean_proximity       15506 non-null  category
 10  people_per_household  15506 non-null  float64 
 11  rooms_per_household   15506 non-null  float64 
 12  bedrooms_per_people   15506 non-null  float64 
 13  bedrooms_per_rooms    15506 non-null  float64 
dtypes: category(1), float64(8), int64(5)
memory usage: 1.7 MB

# Fit Model
lm2 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit()
lm2.summary()

OLS Regression Results
Dep. Variable:	median_house_value	R-squared:	0.743
Model:	OLS	Adj. R-squared:	0.743
Method:	Least Squares	F-statistic:	8956.
Date:	Tue, 18 Jan 2022	Prob (F-statistic):	0.00
Time:	01:09:17	Log-Likelihood:	-1.9073e+05
No. Observations:	15506	AIC:	3.815e+05
Df Residuals:	15500	BIC:	3.815e+05
Df Model:	5
Covariance Type:	nonrobust

	coef	std err	t	P>\|t\|	[0.025	0.975]
Intercept	-5.584e+04	4076.067	-13.700	0.000	-6.38e+04	-4.79e+04
ocean_proximity[T.NEAR COAST]	5.731e+04	1067.209	53.703	0.000	5.52e+04	5.94e+04
median_income	4.88e+04	370.848	131.580	0.000	4.81e+04	4.95e+04
bedrooms_per_rooms	3.69e+05	1.19e+04	30.882	0.000	3.46e+05	3.92e+05
people_per_household	-2.592e+04	562.240	-46.097	0.000	-2.7e+04	-2.48e+04
housing_median_age	924.7246	36.788	25.137	0.000	852.616	996.833

Omnibus:	1156.804	Durbin-Watson:	1.988
Prob(Omnibus):	0.000	Jarque-Bera (JB):	1572.169
Skew:	0.645	Prob(JB):	0.00
Kurtosis:	3.876	Cond. No.	911.

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.

Die Modellperformance liegt nun bei 74,3% und ist somit bisher der beste erzielte Wert der Kennzahl R-squared.
F-statistics hat sich ebenfalls verbessert, auch im Vergleich zu Model 1.
Alle anderen Werte liegen in einem ähnlichen Wertebereich.

fig = sm.graphics.plot_partregress_grid(lm2)
fig.tight_layout(pad=0.2)

fig = sm.graphics.plot_regress_exog(lm2, "housing_median_age")
fig.tight_layout(pad=1.0)

Das Partial Regression Plot und CCPR Plot zeigen eine geringe lineare Korrelation von housing_median_age in Verbindung mit den anderen unabhängigen Variablen auf median_house_value.
Das Residual Plot zeigt eine deutliche Homoskedastizität.

fig = sm.graphics.plot_regress_exog(lm2, "ocean_proximity[T.NEAR COAST]")
fig.tight_layout(pad=1.0)

Die Plots von ocean_proximity[T.NEAR COAST] zeigen keine große Veränderung zu Modell lm.

In der linearen Regression muss eine Colinearität der unabhängigen Variablen unbedingt vermieden werden. Mithilfe des Variance Inflation Factors (VIF) kann der Grad der Colinearität der im Modell verwendeten Variablen berechnet werden. Der kleinstmögliche Wert ist 1. Problematisch wird es ab Werten > 5.

from statsmodels.stats.outliers_influence import variance_inflation_factor
from statsmodels.tools.tools import add_constant

# choose features and add constant
features = add_constant(train_dataset2[['median_income', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']])
# create empty DataFrame
vif = pd.DataFrame()
# calculate vif
vif["VIF Factor"] = [variance_inflation_factor(features.values, i) for i in range(features.shape[1])]
# add feature names
vif["Feature"] = features.columns

vif.round(2)

C:\ProgramData\Anaconda3\lib\site-packages\statsmodels\tsa\tsatools.py:142: FutureWarning:

In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only

	VIF Factor	Feature
0	82.06	const
1	1.80	median_income
2	1.79	bedrooms_per_rooms
3	1.01	people_per_household
4	1.03	housing_median_age

Der Wert der Konstanten kann ignoriert werden. Ansonsten ist die Colinearität der Variablen eher gering.

Die Modelle lm und lm2 sind sehr ähnlich. R-squared ist bei lm2 ca. 1% besser, dafür hat sich der Wert von F-statistics von lm zu lm2 um ca. 1000 verschlechtert. AIC und BIC haben sich nicht wesentlich verändert. Die Werte sind von lm zu lm2 sogar geringfügig kleiner (also besser) geworden, obwohl für zusätzliche Features ein Penalty Term einberechnet wird, um overfitting zu vermeiden. Das Partial Regression Plot der Variablen median_housing_age zeigt zwar nur eine geringe Linearität, dennoch scheint die Aufnahme der Variablen das Modell etwas zu verbessern.

Model 2 - Lasso Regression¶

Nachfolgend wird eine Lasso Regression mit statsmodels fit_regularized() Methode trainiert. Der Wert von alpha wird auf 5 gesetzt.

lasso1 = smf.ols(formula='median_house_value ~ median_income + ocean_proximity + bedrooms_per_rooms + people_per_household + housing_median_age', data=train_dataset2).fit_regularized(method='sqrt_lasso', alpha=5.0, L1_wt=1.0, max_iter=1000)

print("Parameters: ", lasso1.params)

Parameters:  Intercept                        -50616.495719
ocean_proximity[T.NEAR COAST]     57809.976303
median_income                     48394.058880
bedrooms_per_rooms               351976.704469
people_per_household             -26004.829908
housing_median_age                  918.971072
dtype: float64

Da die Summary() Methode für statsmodels fit_regularized() Medelle nicht verfügbar ist, müssen wichtige Kennzahlen selbst berechnet werden.

# calculate the mean 
train_dataset2['average'] = train_dataset2["median_house_value"].mean()

# calculate error (observation - average) and assign it to dataframe
train_dataset2 = train_dataset2.assign(error = (train_dataset2['median_house_value'] - train_dataset2['average']))

# calculate squared error and assign it to dataframe
train_dataset2 = train_dataset2.assign(error_sq = (train_dataset2['median_house_value'] - train_dataset2['average'])**2)


# obtain predicted value from model lasso
train_dataset2['pred'] = lasso1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2'] = train_dataset2['median_house_value'] - train_dataset2['pred']

# square the residuals 
train_dataset2['error_sq_2'] = train_dataset2['error_2']**2

# show df
train_dataset2.head(5)

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms	average	error	error_sq	pred	error_2	error_sq_2
14185	-117.08	32.70	37	2176	418	1301	375	2.8750	98900.0	NEAR COAST	3.469333	5.802667	0.321291	0.192096	196674.726944	-97774.726944	9.559897e+09	157722.078408	-58822.078408	3.460037e+09
6125	-117.91	34.11	20	3158	684	2396	713	3.5250	153000.0	NEAR COAST	3.360449	4.429173	0.285476	0.216593	196674.726944	-43674.726944	1.907482e+09	185009.672924	-32009.672924	1.024619e+09
14095	-117.10	32.75	11	2393	726	1905	711	1.3448	91300.0	NEAR COAST	2.679325	3.365682	0.381102	0.303385	196674.726944	-105374.726944	1.110383e+10	119491.512237	-28191.512237	7.947614e+08
14359	-117.22	32.74	52	1260	202	555	209	7.2758	345200.0	NEAR COAST	2.655502	6.028708	0.363964	0.160317	196674.726944	148525.273056	2.205976e+10	394457.593271	-49257.593271	2.426310e+09
18004	-121.99	37.29	32	2930	481	1336	481	6.4631	344100.0	NEAR COAST	2.777547	6.091476	0.360030	0.164164	196674.726944	147425.273056	2.173421e+10	334928.406586	9171.593414	8.411813e+07

# calculate sum of squared error (which is in case of the mean the total error)
TSS = train_dataset2.error_sq.sum()
# print output
print('Sum of squared error (TSS) of model 1:', TSS)

# Sum of squared residuals (SS_R)
SSR = train_dataset2['error_sq_2'].sum()
print('Sum of squared residuals (SSR) of model 1:',SSR)

Sum of squared error (TSS) of model 1: 170542426261911.88
Sum of squared residuals (SSR) of model 1: 43859574696952.75

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM = TSS - SSR
print('SSM of Lasso Model: ', SSM)

# R_Squared: explained sum of squared residuals
R_squared = SSM / TSS
print('R squared Lasso Model:', R_squared)

SSM of Lasso Model:  126682851564959.12
R squared Lasso Model: 0.7428230871443388

R-squared ist durch die Regularisierung nochmal etwas besser geworden.
Für einen optimalen Wert von alpha müsste eigentlich eine Cross-Validation durchgeführt werden. Darauf wird an dieser Stelle jedoch verzichtet.

Model 3 - Cubic Spline¶

Nun wird eine Cubic Spline unter Verwendung von Patsy und Statsmodels trainiert. Eine Verwendung von mehreren unabhängigen Variablen ist dabei nicht ohne Weiteres möglich. Daher wird nachfolgend lediglich die Variable mit der höchsten Korrelation zu median_house_value verwendet. Dabei handelt es sich um median_income.

X_train = train_dataset2['median_income']
y_train = train_dataset2['median_house_value']
X_test = test_dataset['median_income']
y_test = test_dataset['median_house_value']

from patsy import dmatrix

# Generating cubic spline with 3 knots at 25, 40 and 60
transformed_x = dmatrix(
            "bs(train, knots=(3,6,9), degree=3, include_intercept=False)", 
                {"train": X_train},return_type='dataframe')

# Fitting generalised linear model on transformed dataset
spline1 = sm.GLM(y_train, transformed_x).fit()

spline1.predict()
spline1.summary()

Generalized Linear Model Regression Results
Dep. Variable:	median_house_value	No. Observations:	15506
Model:	GLM	Df Residuals:	15499
Model Family:	Gaussian	Df Model:	6
Link Function:	identity	Scale:	4.6148e+09
Method:	IRLS	Log-Likelihood:	-1.9452e+05
Date:	Tue, 18 Jan 2022	Deviance:	7.1525e+13
Time:	00:43:42	Pearson chi2:	7.15e+13
No. Iterations:	3
Covariance Type:	nonrobust

	coef	std err	z	P>\|z\|	[0.025	0.975]
Intercept	1.081e+05	8237.478	13.122	0.000	9.19e+04	1.24e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[0]	-3.748e+04	1.09e+04	-3.448	0.001	-5.88e+04	-1.62e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[1]	7.416e+04	7398.863	10.023	0.000	5.97e+04	8.87e+04
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[2]	1.575e+05	9980.819	15.782	0.000	1.38e+05	1.77e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[3]	4.47e+05	1.06e+04	41.988	0.000	4.26e+05	4.68e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[4]	3.524e+05	1.97e+04	17.854	0.000	3.14e+05	3.91e+05
bs(train, knots=(3, 6, 9), degree=3, include_intercept=False)[5]	4.232e+05	4.12e+04	10.269	0.000	3.42e+05	5.04e+05

# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline1.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2

# show df
train_dataset2.head(5)

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	...	bedrooms_per_rooms	average	error	error_sq	pred	error_2	error_sq_2	pred_spline	error_2_spline	error_sq_2_spline
14185	-117.08	32.70	37	2176	418	1301	375	2.8750	98900.0	NEAR COAST	...	0.192096	196674.726944	-97774.726944	9.559897e+09	157722.078408	-58822.078408	3.460037e+09	154283.626628	-55383.626628	3.067346e+09
6125	-117.91	34.11	20	3158	684	2396	713	3.5250	153000.0	NEAR COAST	...	0.216593	196674.726944	-43674.726944	1.907482e+09	185009.672924	-32009.672924	1.024619e+09	183299.913013	-30299.913013	9.180847e+08
14095	-117.10	32.75	11	2393	726	1905	711	1.3448	91300.0	NEAR COAST	...	0.303385	196674.726944	-105374.726944	1.110383e+10	119491.512237	-28191.512237	7.947614e+08	96456.816818	-5156.816818	2.659276e+07
14359	-117.22	32.74	52	1260	202	555	209	7.2758	345200.0	NEAR COAST	...	0.160317	196674.726944	148525.273056	2.205976e+10	394457.593271	-49257.593271	2.426310e+09	372578.112664	-27378.112664	7.495611e+08
18004	-121.99	37.29	32	2930	481	1336	481	6.4631	344100.0	NEAR COAST	...	0.164164	196674.726944	147425.273056	2.173421e+10	334928.406586	9171.593414	8.411813e+07	319906.620313	24193.379687	5.853196e+08

5 rows × 23 columns

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Cubic Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Cubic Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Cubic Spline:', R_squared_spline)

Sum of squared residuals (SSR) of Cubic Spline: 71524714796268.75
SSM of Cubic Spline:  99017711465643.12
R squared of Cubic Spline: 0.5806045664764725

print(skm.r2_score(train_dataset2['median_house_value'], train_dataset2['pred_spline']))

0.5806045664764725

R squared ist schlechter geworden als in den vorherigen Modellen. Allerdings liegt dies darin begründet, dass nur ein Feature verwendet wurde. Ein direkter Vergleich ist also nicht möglich.

# Create observations
xp = np.linspace(X_test.min(),X_test.max(), 100)
# Make some predictions
pred = spline1.predict(dmatrix("bs(xp, knots=(3, 6, 9), include_intercept=False)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)

plt.plot(xp, pred, label='Cubic spline with degree=3 (3 knots)', color='orange')
plt.legend();

Model 4 - Natural Spline¶

transformed_x3 = dmatrix("cr(train,df = 3)", {"train": X_train}, return_type='dataframe')

spline2 = sm.GLM(y_train, transformed_x3).fit()

# obtain predicted value from model spline
train_dataset2['pred_spline'] = spline2.predict()

# obtain the residuals from statsmodel (resid)
train_dataset2['error_2_spline'] = train_dataset2['median_house_value'] - train_dataset2['pred_spline']

# square the residuals 
train_dataset2['error_sq_2_spline'] = train_dataset2['error_2_spline']**2

# show df
train_dataset2.head(5)

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	...	bedrooms_per_rooms	average	error	error_sq	pred	error_2	error_sq_2	pred_spline	error_2_spline	error_sq_2_spline
14185	-117.08	32.70	37	2176	418	1301	375	2.8750	98900.0	NEAR COAST	...	0.192096	195254.121813	-96354.121813	9.284117e+09	157516.967720	-58616.967720	3.435949e+09	153409.972829	-54509.972829	2.971337e+09
6125	-117.91	34.11	20	3158	684	2396	713	3.5250	153000.0	<1H OCEAN	...	0.216593	195254.121813	-42254.121813	1.785411e+09	182232.226163	-29232.226163	8.545230e+08	180983.759544	-27983.759544	7.830908e+08
14095	-117.10	32.75	11	2393	726	1905	711	1.3448	91300.0	NEAR COAST	...	0.303385	195254.121813	-103954.121813	1.080646e+10	122387.245306	-31087.245306	9.664168e+08	90298.835046	1001.164954	1.002331e+06
14359	-117.22	32.74	52	1260	202	555	209	7.2758	345200.0	NEAR COAST	...	0.160317	195254.121813	149945.878187	2.248377e+10	394045.429765	-48845.429765	2.385876e+09	352306.649672	-7106.649672	5.050447e+07
18004	-121.99	37.29	32	2930	481	1336	481	6.4631	344100.0	<1H OCEAN	...	0.164164	195254.121813	148845.878187	2.215510e+10	331461.374605	12638.625395	1.597349e+08	313580.198527	30519.801473	9.314583e+08

5 rows × 23 columns

# Sum of squared residuals (SSR)
SSR_spline = train_dataset2['error_sq_2_spline'].sum()
print('Sum of squared residuals (SSR) of Natural Spline:',SSR_spline)

# Explained sum of squares  (SS_M = TSS - SS_R)
SSM_spline = TSS - SSR_spline
print('SSM of Natural Spline: ', SSM_spline)

# R_Squared: explained sum of squared residuals
R_squared_spline = SSM_spline / TSS
print('R squared of Natural Spline:', R_squared_spline)

Sum of squared residuals (SSR) of Natural Spline: 67847658465067.86
SSM of Natural Spline:  95453782865582.95
R squared of Natural Spline: 0.5845250482040099

R-squared ist im Vergleich zur Cubic Spline nochmal gesunken. Auch hier wird im Modell nur eine unabhängige Variable verwendet.

# Make predictions
pred = spline2.predict(dmatrix("cr(xp, df=3)", {"xp": xp}, return_type='dataframe'))

# plot
sns.scatterplot(x=X_train, y=y_train)
plt.plot(xp, pred, color='orange', label='Natural spline with df=3')
plt.legend()

<matplotlib.legend.Legend at 0x1a182d14490>

Scikit Learn¶

Data Preprocessing Pipeline¶

Initialisierung der Preprocessing Pipeline. Mithilfe einer Preprocessing Pipelines müssen die Daten nicht manuell transformiert werden, da diesse Aufgaben von der Pipeline übernommen werden.
Für numerische Variablen gehört dazu die Normalisierung der Variablen, sowie die Belegung nicht befüllter Datensätze mit dem Durchschnittswert des Features.
Bei kategorialen Variablen werden fehlnde Werte durch eine Konstante belegt und automatisiert Dummy Variablen für die verschiedenen Ausprägungen der Variablen erzeugt.

# for numeric features
numeric_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='median')),
    ('scaler', StandardScaler()) #Normalisierung (X=(X-Mittelwert) / Standardabweichung)
    ])

# for categorical features  
categorical_transformer = Pipeline(steps=[
    ('imputer', SimpleImputer(strategy='constant', fill_value='missing')), #Konstante bei fehlenden Werten reinschreiben
    ('onehot', OneHotEncoder(handle_unknown='ignore')) #pro Featureausprägung eine Spalte, zutreffendes Feature hat dann Wert 1 und die anderen Wert 0
    ])

# Pipeline
preprocessor = ColumnTransformer(transformers=[
    ('num', numeric_transformer, selector(dtype_exclude="category")),
    ('cat', categorical_transformer, selector(dtype_include="category"))
        ])

Datenvorbereitung¶

Für Scikit Learn muss die abhängige Variable von den unabhängigen getrennt werden. Als Datenbasis wird hier auf train_dataset2 zugegriffen. In diesem Datenset wurden die Outlier nach Cook’s Distance bereits entfernt. Als Features sollen die durch die Filter Methode identifizierten und bereits in statsmodels verwendeten Features genutzt werden.

# create label y
y = train_dataset2['median_house_value']

# create features
X= train_dataset2[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Scikit Learn bietet eine eigene Methode zum Datensplitting. Das Datenset train_dataset2 aus statsmodels wird nun nochmal in Test- und Trainingsdaten aufgeteilt. Dadurch ist Validierung der Modelle bereits mit diesem Testset möglich, bevor am Ende dieses Notebooks mit den bereits zu Beginn zur Seite gelegten Testdaten das beste Modell evaluiert wird.

X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=10)

print("Training size:", len(X_train))
print("Test size:", len(X_test))

Training size: 12404
Test size: 3102

Model 5 - Simple Regression¶

Nun wird ein Pipeline Object erzeugt, welches zunächst die Data Preprocessing Pipeline ausführt und anschließend das Regressionsmodell.

# Create pipeline with model
lm_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lm', LinearRegression())
                        ])

# Fit model
lm_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('lm', LinearRegression())])

In der Variable features_names werden die Namen der Features gespeichert.

feature_names = np.concatenate((X_train.drop(['ocean_proximity'], axis=1).columns.to_numpy(), lm_pipe.named_steps['preprocessor'].transformers_[1][1]['onehot'].get_feature_names_out()))
feature_names

array(['median_income', 'bedrooms_per_rooms', 'people_per_household',
       'housing_median_age', 'x0_INLAND', 'x0_NEAR COAST'], dtype=object)

Der eben erzeugte feature_names Array kann nun mit den Koeffizienten des Modells ergänzt werden, um die Koeffizienten je Feature auszugeben.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lm_pipe.named_steps['lm'].coef_, feature_names)))

Koeffizienten und Namen: 
 [(85229.2089681146, 'median_income'), (19379.143797351342, 'bedrooms_per_rooms'), (-19642.68963026445, 'people_per_household'), (11314.894820139021, 'housing_median_age'), (-28646.670935754682, 'x0_INLAND'), (28646.670935754682, 'x0_NEAR COAST')]

Das trainierte Modell wird nun für die Berechnung der median_house_values der Trainingsdaten verwendet.

y_pred = lm_pipe.predict(X_train)
y_pred

array([186053.02255921,  51495.90700283, 263988.77300622, ...,
       152104.43702161, 122857.34230185, 120051.11769621])

# Metrics trainings set
print('R squared training set:', round(skm.r2_score(y_train, y_pred)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, y_pred), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, y_pred, squared=False), 2))

R squared training set: 74.18
MSE training set: 2816155847.98
RMSE training set: 53067.47

Das Modell hat einen ähnlichen R-squared wie die OLS Regression in statsmodels. Der Root Mean Squared Error des Modells beträgt 50968,88$. Das bedeutet, dass der eigentliche Wert von median_house_value im Schnitt um diesen Betrag von dem durch das Modell berechneten Ergebnis abweicht.

Im Folgenden wird das Modell mit den im Scikit Learn Datensplit erzeugten Validierungsdaten getestet.

# Validation data
pred = lm_pipe.predict(X_test)
r2_test = round(skm.r2_score(y_test, pred)*100, 2)
print('R squared validation set:', r2_test)
mse_test =skm.mean_squared_error(y_test, pred)
print('MSE validation set', round(mse_test, 2))
print('RMSE training set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 74.67
MSE validation set 2877367258.45
RMSE training set: 53641.1

Der R-squared Wert hat sich mit den Validierungsdaten verbessert. Dies ist sehr ungewöhnliche, da unter Verwendung unbekannter Variablen die Metriken normalerweise etwas schlechtere Werte einnehmen. Dies könnte allerdings darin begründet sein, dass mit diesen Validierungsdaten bereits EDA, Feature Engineering und auch Cook’s Distance Outlier Entfernung durchgeführt wurde.
Der durchschnittliche Fehler MSE und auch die durchschnittliche Abweichung von berechneten zu tatsächlichen Modellen sind etwas schlechter geworden. Allerdings nur geringfügig. Das Modell wird daher als gut erachtet.

# get absolute values of coefficients
importance = np.abs(lm_pipe.named_steps['lm'].coef_)

sns.barplot(x=importance, 
            y=feature_names);

Model 6 - Lasso¶

vorgegebenes Alpha¶

In diesem Unterkapitel wird eine Lasso Regression für das Modell durchgeführt. Der Wert für alpha wird initial mit 1 belegt.

# Create pipeline with model
lasso_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso', Lasso(alpha=1, max_iter=10000))
                        ])

lasso_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('lasso', Lasso(alpha=1, max_iter=10000))])

# Metrics trainings set
pred_train = lasso_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))

R squared training set: 74.18
MSE training set: 2816155856.76
RMSE training set: 53067.47

R-squared der Lasso Regression hat sich im Vergleich zur Simple Regression im Modell zuvor nicht verändert. Auch MSE und RMSE sind gleich geblieben. Dies ist auf den initial gewählten, sehr kleinen Wert für alpha von 1 zurückzuführen. Ein geringer Wert in Alpha verändert die Höhe der Koeffizienten im Modell kaum. Daher bleiben auch die Berechnungen des Modells unverändert und somit die Metriken.

# Metrics validation set
pred = lasso_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 74.67
MSE validation set: 2877369775.47
RMSE validation set: 53641.12

Auch die Metriken des Validation Set haben sich nicht verändert im Vergleich zur Simple Regression.

# Obtain model coefficients and names
print('Koeffizienten und Namen: \n', list(zip(lasso_pipe.named_steps['lasso'].coef_, feature_names)))

Koeffizienten und Namen: 
 [(85226.87228779327, 'median_income'), (19376.806990272275, 'bedrooms_per_rooms'), (-19641.750378420238, 'people_per_household'), (11314.128700828936, 'housing_median_age'), (-57290.99849238739, 'x0_INLAND'), (2.77139978999994e-10, 'x0_NEAR COAST')]

# get absolute values of coefficients
importance = np.abs(lasso_pipe.named_steps['lasso'].coef_)

sns.barplot(x=importance, 
            y=feature_names);

Das Diagramm visualiert die Feature Importance. Je länger der Balken, desto eine höhere Gewichtung bekommt das Feature im Modell. Die Feature Importance ergibt sich auch den Koeffizienten der Features. median_income ist das wichtigste Feature im Modell, wie an dem großen Balken sowie dem Koeffizienten des Features zu erkennen.
Die Ausprägung NEAR COAST der Variablen ocean_proximity wurde von der Lasso Regression aus der Regression entfernt. Dieses Verhalten ist bei Lasso gewünscht. Durch die L1-Regularisierung mit alpha führt der Algorithmus automatisch eine Feature Selection durch. Dabei werden relevante Features durch höhere Koeffizienten stärker gewichtet und unwichtige Features durch kleinere Koeffizienten schäwcher gewichtet bzw. sogar komplett aus dem Modell entfernt, indem der Koeffizient auf den Wert 0 reduziert wird. Durch diese Methode soll ein overfitting des Modells durch zu viele Variablen verhindert werden.

k-fold cross validation zur Ermittlung von alpha¶

from sklearn.linear_model import LassoCV

# Lasso with 5 fold cross-validation
model_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lassoCV', LassoCV(cv=5, random_state=0, max_iter=10000))
                        ])
# Fit model
model_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('lassoCV', LassoCV(cv=5, max_iter=10000, random_state=0))])

#Show best value of penalization chosen by cross validation:
alpha = model_pipe.named_steps['lassoCV'].alpha_
print('Best alpha chosen by cross validation:', alpha)

Best alpha chosen by cross validation: 78.94891231105018

Auf Basis von k-fold cross validation wurde der beste Wert für Alpha ermittelt. Dieser liegt bei 79,59. Bei k-fold cross validation wird das Trainingsdatenset nochmal in k kleinere Datensets unterteilt. Eines Dieser kleineren Datensets wird zur Seite gelegt und das Modell mithilfe der anderen und einem Wert für Alpha trainiert. Anschließend wird mit dem zur Seite gelegten kleinen Datenset eine Validierung des Modells mit den unbekannten Daten durchgeführt.
Dieser Vorgang wird für jedes der k - Datensets für verschiedene Werte von Alpha wiederholt. Das Alpha mit dem besten Ergebnis bei der Validierung wird hier ausgegeben.

Mit dem zuvor ermittelten, besten Wert von Alpha wird nun noch ein finales Lasso Modell auf Basis aller Trainingsdaten trainiert.

# Set best alpha
lasso_best_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('lasso_best', Lasso(alpha=alpha))
                        ])
lasso_best_pipe.fit(X_train, y_train)
lasso_best = lasso_best_pipe.named_steps['lasso_best']

print('Model coefficients and names: \n', list(zip(lasso_best.coef_, feature_names)))

Model coefficients and names: 
 [(85045.084637978, 'median_income'), (19194.892651470986, 'bedrooms_per_rooms'), (-19568.52049548915, 'people_per_household'), (11254.486632101924, 'housing_median_age'), (-57108.04741004523, 'x0_INLAND'), (1.98856656715685e-10, 'x0_NEAR COAST')]

Durch den größeren Wert von Alpha wurden die Koeffizienten der einzelnen Features etwas verkleinert. median_income hatte zuvor beispielsweise einen Wert von 85260,00. Nun liegt der Wert bei 85079,34.

# get absolute values of coefficients
importance = np.abs(lasso_best.coef_)

sns.barplot(x=importance, 
            y=feature_names);

An der Reihenfolge der Feature Importance hat sich durch die Veränderung von Alpha wenig verändert.

# Metrics trainings set
pred_train = lasso_best_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))

R squared training set: 74.18
MSE training set: 2816210527.95
RMSE training set: 53067.98

# Metrics tests set
pred = lasso_best_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 74.67
MSE validation set: 2877619408.17
RMSE validation set: 53643.45

Da sich die Koeffizienten durch den neuen Wert von Alpha kaum verändert haben, sind auch die Kennzahlen zur Performance Messung des Modells unverändert.

Model 7 - Natural Splines¶

In diesem Modell wird eine scikit learn spline in Kombination mit einer Data Preprocessing Pipeline trainiert.

# use a spline wit 4 knots and 3 degrees with a ridge regressions
spline = make_pipeline(SplineTransformer(n_knots=4, degree=3), 
                       Ridge(alpha=1))

# Integrate Preprocessor Pipeline
spline_pipe = Pipeline(steps=[
    ('preprocessor', preprocessor),
    ('spline', spline)
                        ])
                     
spline_pipe.fit(X_train, y_train)

Pipeline(steps=[('preprocessor',
                 ColumnTransformer(transformers=[('num',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(strategy='median')),
                                                                  ('scaler',
                                                                   StandardScaler())]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E691C0>),
                                                 ('cat',
                                                  Pipeline(steps=[('imputer',
                                                                   SimpleImputer(fill_value='missing',
                                                                                 strategy='constant')),
                                                                  ('onehot',
                                                                   OneHotEncoder(handle_unknown='ignore'))]),
                                                  <sklearn.compose._column_transformer.make_column_selector object at 0x0000019AF0E69700>)])),
                ('spline',
                 Pipeline(steps=[('splinetransformer',
                                  SplineTransformer(n_knots=4)),
                                 ('ridge', Ridge(alpha=1))]))])

# Metrics trainings set
pred_train = spline_pipe.predict(X_train)
print('R squared training set:', round(skm.r2_score(y_train, pred_train)*100, 2))
print('MSE training set:', round(skm.mean_squared_error(y_train, pred_train), 2))
print('RMSE training set:', round(skm.mean_squared_error(y_train, pred_train, squared=False), 2))

R squared training set: 75.09
MSE training set: 2716528957.34
RMSE training set: 52120.33

# Metrics tests set
pred = spline_pipe.predict(X_test)
print('R squared validation set:', round(skm.r2_score(y_test, pred)*100, 2))
print('MSE validation set:', round(skm.mean_squared_error(y_test, pred), 2))
print('RMSE validation set:', round(skm.mean_squared_error(y_test, pred, squared=False), 2))

R squared validation set: 75.52
MSE validation set: 2781732015.02
RMSE validation set: 52742.13

Durch die Verwendung der Spline hat sich das Modell in Bezug auf R-sqared und auch RMSE nochmal verbessert. Das Modell deckt nun knapp 76% der Varianz ab und die durchschnittliche Abweichung der tatsächlichen Werten zu den berechneten Werten beträgt knapp 51000$.

# Create observations
xp = pd.DataFrame(data=np.linspace(X_test.min(),X_test.max(), 100), columns = X_train.drop(['ocean_proximity'], axis=1).columns)
xp = xp.assign(ocean_proximity=lambda xp: 'INLAND')


# Make some predictions
pred_new = spline_pipe.predict(xp)

# plot
sns.scatterplot(x=X_train['median_income'], y=y_train)

plt.plot(xp['median_income'], pred_new, label='Cubic spline with degree=3', color='orange')
plt.legend();

C:\Users\Admin\AppData\Local\Temp/ipykernel_9496/3121680389.py:2: FutureWarning:

Dropping of nuisance columns in DataFrame reductions (with 'numeric_only=None') is deprecated; in a future version this will raise TypeError.  Select only valid columns before calling the reduction.

Evaluierung der besten 2 Modelle mit Test Datenset¶

In diesem Kapitel wird eine Evaluierung der zwei besten Modellen mithilfe der zu Beginn beiseite gelegten Testdaten durchgeführt.
Zunächst müssen die Testdaten transformiert werden, sodass sie in gleicher Form wie die Trainingsdaten vorliegen und von den Modellen verwendet werden können.

# write out a dict with the mapping of old to new
remap_cat_dict = {
    'NEAR OCEAN': 'NEAR COAST',
    'NEAR BAY': 'NEAR COAST',
    'INLAND': 'INLAND',
    '<1H OCEAN': 'NEAR COAST',
    'ISLAND' : 'NEAR COAST'}

test_dataset.ocean_proximity = test_dataset.ocean_proximity.map(remap_cat_dict).astype('category')

y_test = test_dataset['median_house_value']

X_test = test_dataset[['median_income', 'ocean_proximity', 'bedrooms_per_rooms', 'people_per_household', 'housing_median_age']]

Evaluation Statsmodels OLS lm2¶

# Metrics tests set
X_test = X_test.assign(pred_test1=lambda X_test: lm2.predict(X_test))
r2_test1 = round(skm.r2_score(y_test, X_test.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test1, squared=False), 2))

R squared Statsmodels OLS with Test Set: -2.13
MSE Statsmodels OLS with Test Set: 13381618868.21
RMSE Statsmodels OLS with Test Set: 115678.95

R-squared ist in Model lm2 negativ. Somit generalisiert das Modell schlecht auf unbekannte Testdaten und stellt sogar eine Verschlechterung im Vergleich zur Verwendung des Durchschnitts als Modell dar.
Grund dafür könnte ein Overfitting während dem Modellierungsprozess sein.
Root Mean Squared Error und somit eine durchschnittliche Abweichung von tatsächlichem und berechnetem median_house_value in Höhe von 115678,95$ ist ebenso kein zufriedenstellendes Ergebnis.

Die berechneten Werte sollen in einem Scatterplot den tatsächlichen Werten gegenübergestellt werden. Aufgrund der Datenmenge sollen nur einige Datensätze visualisiert werden. Mithilfe eines Data Splittings werden zufällig einige Datenpunkte ausgewählt.

X_notplot, X_plot, y_notplot, y_plot = train_test_split(X_test, y_test, test_size=0.1, random_state=10)

# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test1'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();

Im Plot ist zu erkennen, dass bei bei einem median_income von 2 einen Datenpunkt gibt, dessen median_house_value negativ geschätzt wurde. Dieser Wert ist unplausibel, da ein Haus keinen negativen Wert hat.
Den Koeffizienten des Modells lm2 ist zu entnehmen, dass nur das Feature people_per_household eine negative Abhängigkeit auf median_house_value hat. Die Ursache für einen negativen Wert wird also durch einen hohen Wert im Feature people_per_household verursacht.

lm2.params

Intercept                        -55841.725014
ocean_proximity[T.NEAR COAST]     57312.739989
median_income                     48796.057520
bedrooms_per_rooms               368968.477487
people_per_household             -25917.607221
housing_median_age                  924.724553
dtype: float64

test_dataset.loc[test_dataset['people_per_household'] > 100]

	longitude	latitude	housing_median_age	total_rooms	total_bedrooms	population	households	median_income	median_house_value	ocean_proximity	people_per_household	rooms_per_household	bedrooms_per_people	bedrooms_per_rooms
13034	-121.15	38.69	52	240	44	6675	29	6.1359	225000.0	INLAND	230.172414	8.275862	0.006592	0.183333

Tatsächlich existiert ein Datensatz mit über 100 people_per_household. Dieser ist vermutlich für den negativen R-squared Wert verantwortlich. Daher wird dieser testweise entfernt und R-squared erneut berechnet.

X_test2 = X_test.drop(13034)
y_test2 = y_test.drop(13034)

# Metrics tests set
X_test2 = X_test2.assign(pred_test1=lambda X_test2: lm2.predict(X_test2))
r2_test1 = round(skm.r2_score(y_test2, X_test2.pred_test1)*100, 2)
print('R squared Statsmodels OLS with Test Set:', r2_test1)
print('MSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1), 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test2, X_test2.pred_test1, squared=False), 2))

R squared Statsmodels OLS with Test Set: 61.37
MSE Statsmodels OLS with Test Set: 5063025615.54
RMSE Statsmodels OLS with Test Set: 71154.94

Wie vermutet, lag der Grund für den negativen Wert in diesem Ausreiser. Jedoch darf auf den Testdaten gemäß der Literatur keine EDA oder Feature Engineering betrieben werden.
Das Modell lm2 ist daher kein gutes Modell, welches ausreichend unbekannte Daten generalisiert.

Evaluation Scikit Learn Spline¶

#Calculate R2
X_test = X_test.assign(pred_test2=lambda X_test: spline_pipe.predict(X_test))
r2_test = round(skm.r2_score(y_test, X_test.pred_test2)*100, 2)
print('R squared test set:', r2_test)
mse_test =skm.mean_squared_error(y_test, X_test.pred_test2)
print('MSE test set', round(mse_test, 2))
print('RMSE Statsmodels OLS with Test Set:', round(skm.mean_squared_error(y_test, X_test.pred_test2, squared=False), 2))

R squared test set: 65.25
MSE test set 4552968088.14
RMSE Statsmodels OLS with Test Set: 67475.69

Dieses Modell berechnet aus dem kompletten Test Datenset einen zufriedenstellenden Wert für R-squared.
Auch RMSE mit einem Wert von 67475,69$ ist in Ordnung.

# plot
sns.scatterplot(x=X_plot['median_income'], y=y_plot, label='actual')

plt.plot(X_plot['median_income'], X_plot['pred_test2'], label='predicted', color='orange', marker='.', linestyle='')
plt.legend();

Den Analysen entsprechend ist das beste Modell ist die Spline von Scikit Learn mit einem R-squared der Testdaten von 65,25%.
Gemäß dem Data-Science Lifecycle würde dieses Modell jetzt deployed werden. Dieser Schritt wird allerdings nicht mehr im Rahmen dieses Projekts durchgeführt.

Applied Statistics Project

Regression Project

Contents

Regression Project¶

Python Setup¶

Datenaufbereitung¶

Daten laden¶

Inspect Data¶

Datentransformation¶

Entfernung der null values¶

Entfernung unwichtiger Features¶

Konvertierung von Datentypen¶

Hinzufügen neuer Features¶

Datensplit¶

Deskriptive Statistik¶

Kennzahlenanalyse numerischer Variablen¶

Analyse kategorialer Variablen¶

Exploratory Data Anaysis¶

erster Überblick¶

Korrelationsanalyse¶

Analyse der kategorialen Variable Ocean_Proximity¶

Analyse numerischer Variablen¶

Modellierung¶

Statsmodel¶

Model 1 - OLS Regression¶

Regression Diagnostics¶

Model 2 - OLS Regression¶

Model 2 - Lasso Regression¶

Model 3 - Cubic Spline¶

Model 4 - Natural Spline¶

Scikit Learn¶

Data Preprocessing Pipeline¶

Datenvorbereitung¶

Model 5 - Simple Regression¶

Model 6 - Lasso¶

vorgegebenes Alpha¶

k-fold cross validation zur Ermittlung von alpha¶

Model 7 - Natural Splines¶

Evaluierung der besten 2 Modelle mit Test Datenset¶

Evaluation Statsmodels OLS lm2¶

Evaluation Scikit Learn Spline¶